Université de Paris Programmation M1MA
1er semestre 2021-2021
Contrôle continu
A rendre pour le lundi 12 octobre 2020
Instructions
— L’Exercice 1 sera traité en Python et l’Exercice 2 en R.
— L’Exercice 1 est à rendre sous la forme d’un fichier Jupyter notebook Nom_Prenom.ipynb.
L’Exercice 2 sera à rendre sous la forme d’un fichier source R Markdown Nom_Prenom.Rmd
et de son export sous la forme d’un fichier Nom_Prenom.html ou Nom_Prenom.pdf
— Bien structurer vos codes avec des commentaires permettant de comprendre votre raisonnement.
Chacun des deux fichiers doit pouvoir être exécuté tel quel.
— Le fichier .ipynb doit être organisé autant que possible sous la forme de fonctions, avec une struc-
turation claire et des noms de variables permettant de comprendre votre raisonnement.
— Les deux fichiers doivent être envoyés à [email protected] au plus tard le lundi 12
octobre 2020.
— Si un nombre K de projets présentent des similitudes trop importantes chacun des projets concernés
aure une note divisée par K.
1 Exercice 1 (Pyhton)
1.1 Préliminaires : Estimation Bayésienne
On s’intéresse à une pièce de monnaie biaisée et on voudrait estimer la probabilité p d’obtenir "pile".
On lance la monnaie n = 9 fois et on obtient les données
y = (y1, . . . , yn) = (1, 0, 1, 1, 1, 0, 1, 0, 1) (1)
où 1 indique l’évènement "pile" et 0 l’évènement "face". On décide d’estimer le paramètre p par une
approche bayésienne. Pour cela, on considère p comme une variable aléatoire pour la quelle on fixe une
loi a priori f(p) et on choisit un modèle f(y|p) pour y conditionnellement à p, c’est à dire pour la
vraisemblance des données. Dans cet exemple, la vraisemblance des données est :
f(y|p) = P(y1, . . . , yn|p) =
(
n
k
)
pk(1− p)n−k
où k = 6 est le nombre de "piles" dans y. Après avoir observé les données, on met à jour la loi de p en
calculant la loi a posteriori à l’aide du Théorème de Bayes :
f(p|y) = f(y|p)f(p)∫
f(y|p)f(p)dp . (2)
Dans cet exercice vous serez amené à calculer f(p|y) à l’aide de la méthode dite de la grid approxima-
tion. Cette méthode consiste en les étapes suivantes :
a. On définit une grille de valeurs dans l’intervalle de définition de p, ici [0, 1]. Par exemple on pourra
prendre m points équidistants 0 ≤ p1 ≤ p2 ≤ . . . ≤ pm ≤ 1.
b. On évalue la densité a priori en chaque point de la grille : f(p1), . . . , f(pm).
c. On calcule la vraisemblance des données en chaque point de la grille : f(y|p1), . . . , f(y|pm).
d. A l’aide de tous les points de la grille on approche l’intégrale au dénominateur de l’équation (2) par
la somme
m∑
i=1
f(y|pi)f(pi)(pi − pi−1),
en utilisant la convention p0 = 0.
e. A l’aide de l’équation (2), on évalue la loi à posteriori dans chaque point de la grille, et on obtient
ainsi une estimation de la densité à posteriori.
1
1.2 Questions
1. Coder une fonction approx_int qui permet d’approcher numériquement∫ b
a
f(x)dx par la somme
m∑
i=1
f(ti)(ti − ti−1)
évaluée sur la grille grid = (t0, t1, . . . , tm), telle que ti = a+ (b− a) im , obtenue en discrétisant [a, b]
en m intervalles de même longueur.
Tester cette fonction en l’appliquant f(x) = sin(x) sur [0, pi]. Comparer le résultat obtenu avec la
vraie valeur de l’intégrale et au résultat de la fonction scipy.integrate.romb.
Bonus : Coder la fonction de sorte qu’elle accepte en argument soit la fonction (f, grid) soit
(f, a, b,M).
2. Dans chacun des cas suivant, tracer côte à côte la loi a priori p 7→ f(p), la vraisemblance p 7→ f(y|p)
et la loi a posteriori p 7→ f(p|y) pour les données y données en (1) et les lois a priori p 7→ f(p)
suivantes.
(a) Lorsque f est la densité d’une loi uniforme U[0,1].
(b) Lorsque f est la densité suivante : f(p) = 2× 1 1
2(c) Lorsque f est la densité suivante : f(p) = 52 × exp(−5|p− 0.5|).
On veillera à donner un titre à chacun des graphes, ainsi que préciser les axes des abscisses et
des ordonnées. On pourra créer trois fonctions (qui acceptent des vecteurs en entrée) prior pour
p 7→ f(p), likelihood pour p 7→ f(y|p) et post pour p 7→ f(p|y).
A partir de la loi a posteriori f(p|y) on peut définir un estimateur de p, on peut prendre comme
estimateur de p le mode de la loi a posteriori 1 :
pˆ = argmaxpf(p|y).
On va étudier cet estimateur dans les questions qui suivent.
3. Coder une fonction gd qui trouve un extrema local d’une fonction sur un intervalle I ⊂ R par la
méthode dite de descente de gradient. Cette approche utilise le fait que le maximum est à chercher
parmi les points d’annulation de la dérivée. On pourra aussi utiliser que la dérivée d’une fonction
f au point x0 s’approche par df(x0) = (f(x0 + h)− f(x0 − h)/(2h), pour un choix de h petit. On
en déduit l’algorithme dit de descente de gradient suivant :
[Initialisation. ] On choisit une valeur initiale x0 ∈ I.
[Itérations. ] On met à jour la valeur xk par xk+1 donné par :
xk+1 = xk + γdf(xk),
pour un paramètre γ à calibrer. On pourra prendre γ = 1/100 par défaut.
[Fin. On se fixe un critère d’arrêt, soit petit et on s’arrête lorsque : |df(xk)| ≤ .
Par précaution, on ajoute aussi un nombre d’itération maximal N = 1000.
Expliquer pourquoi cet algorithme permettra de calculer pˆ.
Coder une fonction qui prend en argument une fonction f , un point x0, les paramètres h et γ, la
résolution et renvoie un point d’annulation de la dérivée de f .
Tester votre fonction gd sur la fonction g(x) = 1−x2 et comparer le résultat obtenu avec la fonction
optimize de Python.
Bonus : Discuter les limitations d’un tel algorithme.
4. Calculer l’estimateur pour les différents priors considérés à la question 2. Pour trouver l’argument
maximal de l’a posteriori le dénominateur joue-t-il un rôle ?
Le choix du prior semble-t-il influencer l’estimation ?
1. D’autres estimateurs sont possibles comme la médiane a posteriori par exemple.
2
5. (Bonus) L’estimateur ainsi obtenu vous semble-t-il consistant ? Etablir et mettre en œuvre un
protocole permettant d’offrir une justification numérique.
Pour cela, à la place des données (1) utilisées jusqu’à présent, on pensera à générer un échantillon
Y ∼ B(n, p) pour différentes valeurs de n et p (par exemple p = 0.75). Que constate-t-on avec le
second prior de la question 2 si p = 1/3 ?
2 Exercice 2 (R)
2.1 Préliminaires : Fractions continues
Soit a = (a0, . . . , an) une séquence finie d’entiers strictement positifs, la fraction continue associée à
a est par définition le nombre rationnel
f = a0 +
1
a1 +
1
a2+
1
a3+
1
a4+...
(3)
Par exemple, la fraction rationnelle associée à (3, 1, 1, 5) est
3 +
1
1 + 1
1+ 15
On dira que a est la réduite de f . Pour tout nombre rationnel f il existe une réduite a telle que f admet
la représentation donnée par l’équation (3). Par exemple 15644 a pour réduite la séquence (3, 1, 1, 5) de
l’exemple ci-dessus. Dans cet exercice vous serez amenés à calculer la réduite d’un nombre rationnel à
partir de l’algorithme de division euclidienne. Par souci de simplicité on considérera des rationnels non
négatifs.
2.2 Questions
6. On rappelle que à deux entiers non négatifs a, b, b > 0, on associe de façon unique deux entiers non
négatifs q, r tels que a = qb + r et r < b. Pour calculer le quotient q et le reste r de la division
euclidienne de a par b on pourra soustraire b à a jusqu’au moment où on obtient une différence
inférieure à b : cette différence donnera r et le nombre de soustractions faites donnera q. Plus
précisément, on considère l’algorithme de division euclidienne suivant :
[Entrée. ] Deux entiers a, b avec a ≥ 0 et b > 0
[Initialisation. ] r ← a et q ← 0
[Itérations. ] Tant que r ≤ b : r ← r − b et q ← q + 1
[Sortie. ] Rendre r, q
Coder une fonction div_euclid qui implémente l’algorithme ci-dessus. Dans la fonction, on ajoutera
une commande qui teste si l’argument b mis en entrée est bien strictement positif et qui renvoie le
message d’erreur "Erreur: b doit être strictement positif" sinon. On ajoutera également
une commande qui teste si l’argument a est bien non négatif et qui renvoie le message d’erreur
"Erreur: a ne peut pas être négatif" sinon.
Comme vérification, calculer le quotient et le reste de la division de 77708431 par 2640858 et
comparer ces résultats à ceux obtenus avec les fonctions standard de R.
7. On considère un nombre rationnel non négatif f0f1 , et on cherche la réduite qui lui est associée. On
peut commencer par observer que par division euclidienne il existe deux entiers a0 et f2 < f1 t.q.
f0 = a0f1 + f2. Si f2 = 0 on a f0f1 = a0 et (a0) est la réduite cherché. Si f2 > 0 on a
f0
f1
= a0 +
1
f1
f2
et on peut procéder et trouver par division euclidienne deux entiers a1 et f3 < f2 tels que, si f3 > 0,
f1
f2
= a1 +
1
f2
f3
,
3
ce qui donne
f0
f1
= a0 +
1
a1 +
1
f2
f3
.
On continue de cette sorte en faisant à chaque itération la division euclidienne de fi−1 par fi pour
obtenir ai−1 et fi+1 jusqu’au moment où fn+1 = 0. Par construction, la séquence (a0, . . . , fn) est
la réduite de f0f1 .
Coder une fonction reduite qui prend en argument deux entiers f0 ≥ 0, f1 > 0 et implémente la
procédure décrite ci-dessus pour rendre en sortie la réduite de f0f1 .
Comme vérification on calculera la réduite de 777084312640858 .
8. Coder une fonction frac_cont qui prend en entrée une séquence d’entiers strictement positifs
a = (a0, . . . , an) et rend en sortie la valeur de la fraction continue associé à a en notation décimale
(Rappel : la notation décimale de 12 est 0.5.)
Comme vérification on montrera que la fraction continue associée à la réduite de 777084312640858 a la même
valeur en notation décimale de 777084312640858 .
9. Bonus : Si (a0, . . . an) est la réduite d’un nombre rationnel f , pour tout 0 ≤ p ≤ n la séquence
(a0, . . . , ap) est dite réduite d’ordre p de f . Les fractions continues associées aux réduites d’ordre
p < n donnent des approximations de f .
Pour tout 0 ≤ p < n calculer l’approximation de f = 777084312640858 donnée par la fraction continue
associée à la réduite d’order p de f et la comparer avec la valeur de f . On montrera à l’aide d’un
graphique approprié que ces approximations sont de plus en plus précises pour des valeurs de p de
plus en plus grandes.
4

欢迎咨询51作业君