La densité de la loi Gamma $\Gamma(a,b)$ est donnée par $$f_{(a,b)}(x) = \frac{b^a}{\Gamma(a)}x^{a-1}e^{-bx},\quad x>0$$ avec des paramètres $a>0$ et $b>0$.
Soit $\mathbf x=(x_1,\dots,x_n)$ une réalisation de $\mathbf X=(X_1,\dots,X_n)$ où les $X_i$ sont i.i.d. de la loi Gamma $\Gamma(a,b)$.
L'objectif de ce TP est l'estimation des paramètres $a$ et $b$ à partir des observations $\mathbf x$.
La fonction de log-vraisemblance dans ce modèle s'écrit $$\ell(a,b) = \sum_{i=1}^n\log f_{(a,b)}(x_i)= na\log(b) -n\log(\Gamma(a)) +(a-1)\sum_{i=1}^n\log(x_i) -b\sum_{i=1}^nx_i$$
L'estimateur du maximum de vraisemblance de $a$ et $b$ est donné par $$(\hat a, \hat b) = \arg\max_{a>0,b>0}\ell(a,b).$$
Pour une représentation graphique d'une fonction $f:\mathcal A\subset\mathbb R^2\mapsto \mathbb R$ on peut tracer un graphique en trois dimensions. Plusieurs solutions existent sous R, et chacune nécessistent des libraries spécifiques, qu'il faut éventuellement installer sur votre machine (avec la commande install.packages()). En revanche, sur les machines à l'UTES vous n'avez pas les droits pour installer des packages.
Nous présentons ici deux solutions pour tracer les valeurs d'une fonction $f(x,y)$. La première crée un graphique interactif, très abouti. Cette solution est à préférer. Malheureusement, elle ne fonctionne (toujours) pas sur les machines à l'UTES. La deuxième solution crée un graphique statique, moins joli que l'autre, mais qui fonctionne à l'UTES.
Pour tracer les valeurs d'une fonction $f(x,y)$, on peut utiliser la fonction plot_ly du package plotly :
où x.grid et y.grid sont des vecteurs contenant des grilles de valeurs pour l'abscisse et l'ordonnée, et z.values est une matrice de taille length(x.grid)$\times$length(y.grid) avec les valeurs correspondantes de $f(x,y)$. L'option type='surface' interpole les points pour une représentation continue de la fonction $f$.
Pour tracer les valeurs d'une fonction $f(x,y)$, on peut aussi utiliser la fonction ggplot et pour cela il faut charger les packages ggplot2 et plyr. L'instruction est la suivante :
où df est un dataframe avec trois colonnes, ici de nom x, y et f_xy. Chaque ligne de df donne les coordonnées d'un point à tracer, à savoir $(x,y,f(x,y))$. Par conséquent, pour tracer tous les points de $f$ pour des grilles de valeurs x.grid et y.grid, il faut créer un dataframe df avec length(x.grid)$*$length(y.grid) lignes.
On peut montrer que l'EMV n'admet pas d'expression explicite (voir TD). En revanche, on peut recourir à la méthode de Newton-Raphson pour l'approcher. Plus précisément, on peut l'appliquer afin de trouver des points critiques de $\ell(a,b)$.
Rappelons que la méthode de Newton-Raphson est un algorithme itératif. Notons $(a^{(t)},b^{(t)})^T$ les valeurs actuelles de paramètres. Une itération de l'algorithme consiste à calculer les nouvelles valeurs des paramètres $(a^{(t+1)},b^{(t+1)})^T$ par la formule $$(a^{(t+1)},b^{(t+1)})^T = (a^{(t)},b^{(t)})^T-[H(a^{(t)},b^{(t)})]^{-1}\nabla\ell(a^{(t)},b^{(t)}),\qquad\qquad(*)$$ où $H$ désigne la matrice Hessienne de $\ell$. Pour le gradient et l'inverse de la Hessienne on obtient (voir TD pour les calculs détaillés) \begin{align*} \nabla\ell(a,b) &= \left(\frac\partial{\partial a} \ell(a,b),\frac\partial{\partial b} \ell(a,b)\right)^T =\left( n\log b- n (\log\Gamma(a))'+\sum_{i=1}^n\log x_i,\quad \frac{na}b-\sum_{i=1}^n x_i \right)^T \end{align*} et \begin{align*} [H(a,b)]^{-1} =\left(\begin{array}{cc} \frac{\partial^2}{\partial^2 a} \ell(a,b)&\frac{\partial^2}{\partial a\partial b} \ell(a,b)\\ \frac{\partial^2}{\partial a\partial b} \ell(a,b)&\frac{\partial^2}{\partial^2 b} \ell(a,b) \end{array}\right)^{-1} &= \frac1{n\left(1-a(\log\Gamma(a))''\right)} \left(\begin{array}{cc} a&b\\ b&b^2(\log\Gamma(a))'' \end{array}\right) \end{align*}
Ecrire une fonction R qui effectue une mise à jour des paramètres. Cette fonction prend en argument les valeurs actuelles des paramètres, $a^{(t)}$ et $b^{(t)}$, ainsi que les données. Elle renvoie les nouvelles valeurs des paramètres $a^{(t+1)}$ et $b^{(t+1)}$ obtenues par $(*)$. Pour les dérivées de $\log\Gamma(a)$ utiliser les fonctions digamma et trigamma (regarder l'aide de R pour les détails!!). On peut compléter la fonction suivante :
La méthode de Newton-Raphson consiste à répéter les mises à jour décrites ci-dessus pour $t=1,2,\dots$ par une boucle while jusqu'à convergence. Il nous faut alors un critère d'arrêt qui vérifie après chaque itération si l'algorithme a convergé.
Plusieurs critères d'arrêt sont envisageable. Ici nous allons considérer le critère qui repose sur la fonction à maximiser $\ell(a,b)$. Ainsi, on arrête les calculs dès que la valeur de la log-vraisemblance est stalbe. Plus précisément, on arrête dès que $$\left|\frac{\ell(a^{(t+1)},b^{(t+1)})-\ell(a^{(t)},b^{(t)})}{\ell(a^{(t+1)},b^{(t+1)})}\right|<\varepsilon,$$ où $\varepsilon>0$ est un seuil fixé.
La méthode de Newton nécessite un point d'initialisation $(a^{(0)},b^{(0)})^T$. C'est à l'utilisateur de le choisir. Nous verrons l'importance d'un bon choix du point initial.
Ecrire une fonction, nommée emv_gamma, qui calcule l'EMV en utilisant la méthode de Newton-Raphson. Elle prend en argument les données et les valeurs initiales $a^{(0)}$ et $b^{(0)}$. En partant de ces valeurs initiales on itère des mises à jour des paramètres selon la méthode de Newton jusqu'à convergence. Utilisez le critère d'arrêt ci-dessus avec le seuil $\varepsilon=10^{-3}$. La fonction renvoie les estimateurs de $a$ et $b$ et le nombre d'itérations effectuées. Vous pouvez compléter la fonction suivante :
Afin d'éviter que l'algorithme tourne indéfiniment, modifier la fonction en sorte que l'on s'arrête lorsque l'algorithme à convergé ou après au plus 100 itérations.
Testez votre fonction sur des données simulées. Testez différentes valeurs de paramètres $a$ et $b$ ainsi que différentes tailles d'échantillons. Pour commencer, initialisez toujours avec les vraies valeurs des paramètres. Dans ce cas, l'algorithme doit converger très rapidement. Si ce n'est pas le cas, votre programme contient des bugs qu'il faut corriger !
Un autre estimateur des paramètres $a$ et $b$ dans ce modèle est donné par l'estimateur par la méthode des moments (EMM). En TD vous montrez que cet estimateur est donné par $$\tilde a = \frac{(\bar X_n)^2}{s_X^2},\qquad \tilde b = \frac{\bar X_n}{s_X^2},$$ où $s_X^2=\frac1n\sum_{i=1}^n(X_i-\bar X_n)^2$ désigne la variance empirique.
Ecrire une fonction nommée emm_gamma qui prend en argument les observations et renvoie l'estimateur par la méthode des moments de $a$ et $b$.
Modifier votre fonction emv_gamma en sorte que par défaut elle utilise comme valeurs initiales l'estimateur par la méthode des moments $\tilde a$ et $\tilde b$. Etudier de nouveau la convergence de l'algorithme sur des données simulées.
Charger le package datasets (par l'instruction library(datasets)) qui contient le jeu de données nommé lynx. Estimer les paramètres de la loi Gamma pour les données des lynx. Superposer la densité estimée à l'histogramme. Comparer avec la densité exponentielle ajustée aux données. Commenter.
On veut comparer les deux estimateurs, EMV et EMM, des paramètres de la loi Gamma. On cherche à savoir lequel a le plus petit risque quadratique.
Rappelons que le risque quadratique d'un estimateur $\hat\theta$ d'un paramètre $\theta$ est défini par $\mathcal R(\hat\theta,\theta)= \mathbb E[\|\hat\theta-\theta\|^2]$.
Ecrire un programme qui
Faites des simulations pour savoir quel estimateur a la meilleur performance.