Estimation de paramètres de la loi Gamma

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$.

Estimateur du maximum de vraisemblance

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).$$

Graphiques en 3D

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.

1. Plot 3D interactif

Pour tracer les valeurs d'une fonction $f(x,y)$, on peut utiliser la fonction plot_ly du package plotly :

plot_ly(x=x.grid, y=y.grid, z=z.values, type="surface")

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$.

2. Plot 3D statique

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 :

ggplot(df, aes(x=x, y=y, z=f_xy)) + geom_contour(bins = 60, aes(colour = ..level..))

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.

Question 1

  • Ecrire une fonction qui évalue la fonction de log-vraisemblance $\ell(a,b)$. Elle prend en argument des valeurs pour les paramètres $a$ et $b$ et les observations $\mathbf x=(x_1,\dots,x_n)$. Pour calculer la fonction Gamma on utilisera la fonction gamma.
  • Ecrire une fonction qui permet de visualiser la fonction de log-vraisemblance en trois dimensions. Créer un plot interactif si votre machine le permet.
  • Faites des représentations graphiques pour différents jeux de données simulés de tailles différentes et avec différentes valeurs de paramètres $a$ et $b$ (la fonction rgamma tire des observations de la loi Gamma). Quelle est l'allure typique de la fonction de log-vraisemblance ? Combien de points critiques a-t-elle ? Ou est situé le maximum de la fonction ?

Calcul de l'EMV par la méthode de Newton-Raphson

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)$.

Une mise à jour de Newton-Raphson

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*}

Question 2

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 :

update_newton <- function(param,obs){ # param - liste de paramètres: param$a, param$b # obs - vecteur de données # calcul des nouvelles valeurs des paramètres : ... a.new <- ... b.new <- ... return(list(a=a.new,b=b.new)) }

Itérer les mises à jour

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é.

Critère d'arrêt

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é.

Initialisation

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.

Question 3

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 :

emv_gamma <- function(obs,param){ # obs - vecteur de données # param - liste de paramètres: param$a, param$b param.old <- param not.converged <- TRUE ... while(not.converged){ # mise à jour des paramètres : param.new <- update_newton(param,obs) # vérifier le critère d'arrêt : not.converged <- ... param.old <- param.new ... } return(list(param=param, nb.iter=nb.iter)) }

Question 4

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.

Question 5

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 !

Question 7

  • Essayez maintenant d'initialiser avec des points initiaux de plus en plus éloignés des vraies valeurs des paramètres. Qu'observez-vous ?
  • On observe parfois que cela produit des erreurs. En fait, il arrive que la fonction update_newton renvoie des valeurs négatives pour a et/ou b (ce qui n'est pas admissible du point de vue d'interprétation de $a$ et $b$ comme paramètres de la loi Gamma, mais la méthode de Newton ne respecte pas le domaine de définition). Dans ce cas, le calcul de la fonction de log-vraisemblance produit des erreurs. Modifiez la fonction emv_gamma en sorte qu'elle arrête l'algorithme dès que update_newton renvoie des valeurs négatives.

Initialisation avec l'estimateur par la méthode des moments

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.

Question 8

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$.

Question 9

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.

Application aux données des lynx

Question 10

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.

Comparaison de l'EMV et de l'EMM

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]$.

Question 11

Ecrire un programme qui

  • génère un grand nombre de jeux de données de loi Gamma,
  • évalue les deux estimateurs sur ces données et
  • calcule une approximation du risque quadratique des deux estimateurs.

Question 12

Faites des simulations pour savoir quel estimateur a la meilleur performance.