Simulations Monte-Carlo

L'objectif de ce TP est l'étude des propriétés d'un estimateur par des simulations.

Vous allez vous habituer à programmer des fonctions R. N'hésitez pas à consulter le mini-poly sur la Programmation en R si besoin.

Modèle statistique

Tout au long de ce TP, nous allons considérer le modèle statistique suivant. Soit $X$ la variable aléatoire définie comme $$X=Z+\delta, \qquad\qquad(*)$$ où $\theta\in\mathbf R$ est une constante et $Z$ une variable aléatoire de loi de Student $t_q$ à $q$ degrés de liberté.

Nous supposons que $q$ est connu et le paramètre $\theta$ est inconnu. Le but est l'estimation de $\theta$ à partir de $n$ copies i.i.d. $(X_1,\dots,X_n)$ de $X$.

Remarquons qu'il s'agit d'un modèle de translation et $\theta$ est un paramètre de position.

Exercice 1

  1. Ecrire une fonction rttrans qui génère un échantillon $(X_1,\dots,X_n)$ selon le modèle $(*)$. Cette fonction prend en argument
    • n : la taille d'échantillon,
    • q : le degré de liberté de la loi de Student de $Z$ et
    • theta : le paramètre de position $\theta$.
      Elle retourne un vecteur de taille n avec un échantillon simulé sous le modèle $(*)$ avec les paramètres q et theta. On pourra utiliser la fonction rt pour générer des réalisations de la loi de Student (centrée).
  2. Vérifier que votre fonction rttrans fonctionne. (Pour cela il faut soit enregistrer et sourcer votre script, soit sélectionner les lignes de votre fonction et les exécuter une fois. Cela crée l'objet rttrans que vous pourrez ensuite utiliser comme toute autre fonction prédéfinie sous R.)

Moyenne empirique

Pour tout $q>1$ la loi de Student $t_q$ est intégrable et donc $X$ l'est aussi vérifiant $\mathbb E[X]=\theta$. Par la méthode des moment, on obtient la moyenne empirique $\bar X_n$ comme estimateur de $\theta$.

On aimerait étudier la qualité de cet estimateur en fonction de la taille d'échantillon $n$, la valeur du paramètre $\theta$ et le degré de liberté $q$.

Risque quadratique

Rappelons qu'un critère pour évaluer la qualité d'un estimateur $\hat \theta$ de $\theta$ est son risque quadratique définit par $$\mathcal R(\hat\theta,\theta) = \mathbb E_{\theta}\left[(\hat\theta-\theta)^2\right].$$

Dans notre modèle le risque quadratique de la moyenne empirique $\bar X_n$ n'est pas explicite, car la loi de la variable aléatoire $T =(\bar X_n-\theta)^2$ est difficile à déterminer. Nous choisissons alors d'analyser le comportement de $\bar X_n$ sur des données simulées.

Plus précisément, l'idée consiste à approcher l'espérance dans la définition du risque quadratique par une moyenne empirique en simulant un grand nombre de réalisations $T_k$ de la variable $T=(\hat\theta-\theta)^2$. Comment faire puisque la loi de $T$ étant inconnue ? En fait, l'unique aléa dans $T$ vient de l'échantillon $(X_1,\dots,X_n)$, car l'estimateur $\hat\theta$ est une fonction mesurable sur l'espace des observations $(X_1,\dots,X_n)$. Ainsi, $T=(\hat\theta-\theta)^2= g(X_1,\dots,X_n)$ pour une fonction mesurable $g$. Par conséquent, on obtient une réalisation $T_k$ de $T$ :

  1. en générant un échantillon $(x^{(k)}_1,\dots,x^{(k)}_n)$ (c'est facile avec rttrans) et
  2. en évaluant $g$ en cet échantillon simulé : $T_k = g(x^{(k)}_1,\dots,x^{(k)}_n).$

On répète cette démarche $K$ fois pour créer un échantillon $(T_1,\dots,T_K)$. Ensuite, il ne reste à calculer la moyenne empirique $\bar T_K$ des $(T_1,\dots,T_K)$ qui est une approximation du risque quadratique.

C'est le principe des simulations dites de Monte Carlo. Ce procédé est justifié par la loi des grands nombres, car $$\bar T_k \stackrel{P}{\longrightarrow} \mathcal R(\hat\theta,\theta),\qquad K\to\infty.$$ Ainsi, plus le nombre $K$ de réalisations $T_k$ est grand, mieux est l'approximation du risque quadratique $\mathcal R(\bar X_n,\theta)$ par $\bar T_K$. On appele $\bar T_k$ le risque quadratique empirique.

Exercice 2

  1. Ecrire une fonction R, nommée moy_student, qui prend en argument :
    • n : la taille d'échantillon,
    • q : le degré de liberté et
    • theta : la valeur du paramètre.
      La fonction moy_student simulera un échantillon $(X_1,\dots,X_n)$ de taille n selon le modèle $(*)$ avec les paramètres q et theta. Ensuite, elle calculera la moyenne empirique $\bar X_n$ qui sera renvoyée en sortie.
  2. Ajouter dans la fonction moy_student l'argument K qui est le nombre d'échantillons à générer et dans le corps de la fonction rajouter une boucle for afin de générer pas un, mais K échantillons de taille n. La fonction renverra un vecteur avec les K moyennes empiriques $\bar X_n^{(k)}$.

Exercice 3

Ecrire une fonction nommée risque qui prend en argument une valeur theta ainsi qu'un vecteur hat.theta avec des estimés de theta. La fonction renverra le risque quadratique empirique associé.

Exercice 4

  1. Fixons n=50, q=5 et theta=0. A partir de quelle valeur de K le risque quadratique de $\bar X_n$ reste relativement stable ? On utilisera cette valeur de K pour toutes les simulations suivantes. (Dans la définition de moy_student on pourra définir cette valeur comme la valeur par défaut de K.)
  2. Maintenant on veut analyser la performance de l'estimateur $\bar X_n$ en fonction de la taille d'échantillon $n$. On fixera les paramètres q=5 et theta=0. Ecrire un script qui
    • appelle la fonction moy_student pour des différentes tailles d'échantillon (on prendra n = 10, 50, 100, 200, 1000),
    • on crée un dataframe pour y stocker les estimées $\bar X_n$ renvoyées par moy_student,
    • on calcule les risques quadratiques associés pour les différentes valeurs de n et
    • on trace les boxplots des estimées $\bar X_n$ pour les différents n.
      Interpréter les résultats.
  3. Afin d'analyser maintenant l'impact de la valeur de theta sur l'estimateur $\bar X_n$, on fixe les paramètres q=5 et n=50 et on applique la même démarche qu'à la question précédente en variant la valeur de theta (on peut considérer theta = 0, 10, 300, -4444). On calculera alors les différents risques quadratiques, et on tracera les boxplots des différences $\bar X_n-\theta$ (au lieu de $\bar X_n$). Interpréter les résultats.
  4. Même question pour le degré de liberté q. Fixons n=50 et theta=0. Calculer le risque quadratique et tracer le boxplot des $\bar X_n$ pour les valeurs suivantes de q : 1, 2, 5, 10, 50. Qu'observez vous ? Expliquez.
In [ ]: