Université de Paris
1er semestre 2021-2021
Programmation M1MA
— —
— — — —
1
1.1
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 à celine.duval@parisdescartes.fr 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.
Exercice 1 (Pyhton)
Préliminaires : Estimation Bayésienne
Contrôle continu
A rendre pour le lundi 12 octobre 2020
Instructions
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 :
n k n−k f(y|p)=P(y1,…,yn|p)= k p (1−p)
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 (y|p)f (p)
f(p|y) = 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
f(y|pi)f(pi)(pi −pi−1), i=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 m f(x)dx par la somme
f(ti)(ti − ti−1)
évaluée sur la grille grid = (t0,t1,…,tm), telle que ti = a+(b−a) i , obtenue en discrétisant [a,b]
m
en m intervalles de même longueur.
Tester cette fonction en l’appliquant f(x) = sin(x) sur [0,π]. 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 → f(p), la vraisemblance p → f(y|p) et la loi a posteriori p → f(p|y) pour les données y données en (1) et les lois a priori p → 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 × 11
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. ]Deuxentiersa,baveca≥0etb>0 [Initialisation. ] r ← a et q ← 0
[Itérations. ]Tantquer≤b:r←r−betq←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 f0 , et on cherche la réduite qui lui est associée. On f1
peut commencer par observer que par division euclidienne il existe deux entiers a0 et f2 < f1 t.q. f0=a0f1+f2.Sif2=0onaf0 =a0et(a0)estlaréduitecherché.Sif2>0ona
f1
f0 = a0 + 1
f1 f1
f2
et on peut procéder et trouver par division euclidienne deux entiers a1 et f3 < f2 tels que, si f3 > 0,
f1 = a1 + 1 ,
a+1
3 a4 +…
f2
f2 f3
3
ce qui donne
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 f0 . f1
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 f0 .
f1
Comme vérification on calculera la réduite de 77708431 . 2640858
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 1 est 0.5.) 2
Comme vérification on montrera que la fraction continue associée à la réduite de 77708431 a la même
valeur en notation décimale de 77708431 . 2640858
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 = 77708431 donnée par la fraction continue 2640858
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.
f0=a0+ 1 . f a+1
1 1f2 f3
2640858
4