Données simulées

Nous allons commencer par simuler des données selon une densité mélange. Soit \(Z_1,...,Z_n\sim_{i.i.d.}Z\) où \(Z\) suit la loi uniforme sur \(\{1,2,3\}\).Puis

\[ X|Z=j\sim\mathcal N(\mu_j,\sigma_j^2),\qquad j=1,2,3 \] où \(\mu_j\), \(j=1,2,3\) et \(\sigma_j^2\), \(j=1,2,3\) sont des paramètres sur lesquels on pourra jouer.

La densité de \(X\) suit une loi mélange de gaussiennes` \[ f_X(x) = \sum_{j=1}^3\frac13f_{\mu_j,\sigma_j^2}(x), \] où \(f_{\mu_j,\sigma_j^2}\) est la densité de la loi \(\mathcal N(\mu_j,\sigma_j^2)\).

  1. Simuler un vecteur \(Z\) de longueur \(n=100\) contenant une réalisation de \(Z_1,...,Z_n\). On pourra utiliser la fonction sample.int() avec l’option replace=TRUE.

  2. Simuler maintenant un vecteur \(X\) de même longueur contenant une réalisation de \(X_1,...,X_n\). On pourra prendre \(\mu_1=-3\), \(\mu_2=0\), \(\mu_3=1\), \(\sigma_1^2=1\), \(\sigma_2^2=1\), \(\sigma_3^2=0.03\).

mu1 = -3
mu2 = 0
mu3 = 1
sigma12 = 1
sigma22 = 1
sigma32 = 0.03

X = rnorm(n,mu1,sqrt(sigma12))
X[Z==2] = rnorm(sum(Z==2),mu2,sqrt(sigma22))
X[Z==3] = rnorm(sum(Z==3),mu3,sqrt(sigma32))
  1. À l’aide de la fonction density(), tracer un premier estimateur à noyau de la densité de \(X\). On pourra ajouter les points de \(X\) sous la courbe à l’aide de la commande rug(X). Tracer également la vraie densité de \(X\) et comparer.

  1. Vérifier dans l’aide de la fonction density() comment est sélectionné la fenêtre. Essayer de modifier la fenêtre pour améliorer l’estimation.

  2. La fonction density a d’autres méthodes de sélection de la fenêtre. Essayer les options bw='SJ', bw='ucv' et bw='bcv'.

6.(bonus) Implémenter la méthode PCO et comparer. On prendra par exemple comme ensemble de fenêtres \[ \mathcal H_n = \left\{\frac{k^2\|K\|_1\|K\|_{\infty}}{n}, k=1,...,k_{\max}\right\} \] avec \(k_{\max}=20\).

Pour calculer la valeur théorique de \(\|K\|_1\) et \(\|K\|_\infty\) on pourra remarquer dans un premier temps que le noyau par défaut est le noyau gaussien \[ K(t)=\frac1{\sqrt{2\pi}}e^{-t^2/2}. \] Pour une fonction \(f\), on approchera la norme \(\|f\|^2_2\) par \[ \frac1{p}\sum_{j=1}^pf^2(t_j), \] pour \(t_1,....,t_p\) tels que \(t_j=(j-1)/(p-1)\) et le produit scalaire \(\langle f,g\rangle\) par \[ \frac1p\sum_{j=1}^pf(t_j)g(t_j). \] Attention à prendre \(p\) grand, par exemple \(p=5000\) pour que la méthode d’approximation par la somme soit précise même lorsque l’estimateur est irrégulier (ce qui est le cas notamment pour les plus petites valeurs de \(h\)).

Données réelles

  1. À l’aide de ce que vous avez appris sur les estimateurs à noyau de la densité. Proposez un estimateur à noyau pour la densité des observations :

On pourra comparer et commenter.