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.
n=100
Z = sample.int(3,n,replace=T)
  1. 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.
plot(density(X,bw=0.1),ylim=c(0,0.85),main='Estimation de la densité')
rug(X)
x = seq(-6,2,length.out=100)
points(x,(dnorm(x,mu1,sqrt(sigma12))+dnorm(x,mu2,sqrt(sigma22))+dnorm(x,mu3,sqrt(sigma32)))/3, type='l',col='darkgreen')
legend('topleft',c('Vraie densité','Estimation'),lty=c(1,1),col=c('darkgreen',1))

# la valeur 0." semble donner une meilleure estimation.
  1. 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'.
plot(density(X,bw='SJ'),ylim=c(0,0.85),main='Estimation de la densité')
lines(density(X,bw='ucv'),col=2)
lines(density(X,bw='bcv'),col=4)
## Warning in bw.bcv(x): minimum occurred at one end of the range
rug(X)
x = seq(-6,2,length.out=100)
points(x,(dnorm(x,mu1,sqrt(sigma12))+dnorm(x,mu2,sqrt(sigma22))+dnorm(x,mu3,sqrt(sigma32)))/3, type='l',col='darkgreen')
legend('topleft',c('Vraie densité','SJ','ucv','bcv'),lty=rep(1,4),col=c('darkgreen',1,2,4))

# la validation croisée non biaisée (UCV) semble donner les meilleurs résultats..

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

kmax = 50
H = (1:kmax)^2/(n*sqrt(2*pi))
hmin = min(H)
p=5000
fhat = matrix(NA,kmax,p)
fhatmin = density(X,bw=hmin,n=p)
x = fhatmin$x
fhat[1,] = fhatmin$y
for (j in 2:kmax){
  fhat[j,] = density(X,bw=H[j],n=p)$y
}

# Calcul du critere PCO
b = rep(NA,kmax)
v = rep(NA,kmax)
crit = rep(NA,kmax)
for (k in 1:kmax){
  b[k] = mean((fhat[k,]-fhat[1,])^2)
  v[k] = 2*mean(dnorm(x,sd = hmin)*dnorm(x,sd = H[k]))/n
  crit[k]  = b[k]+v[k]
}
khat = which.min(crit)
hhat = H[khat]
fhhat = fhat[khat,]
plot(x,fhhat,main = paste("Minimum atteint pour h=",signif(hhat,2)),xlab='t',ylab=expression(hat(f)[hat(h)]),type='l')
points(x,(dnorm(x,mu1,sqrt(sigma12))+dnorm(x,mu2,sqrt(sigma22))+dnorm(x,mu3,sqrt(sigma32)))/3,type='l',col='darkgreen',lwd=2)

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 :
plot(density(precip,bw='ucv'))
rug(precip)

# Les données sont bimodales. On décèle deux groupes de villes : un groupe avec des précipitations situées autour de 15 inches avec une grande dispersion et un autre groupe avec des précipitations moyennes situées autours de 40 inches.
plot(density(faithful[,1],bw='ucv'))
rug(faithful[,1])

# On distingue là aussi deux groupes : un groupe de geysers ayant un temps d'éruption autour de 2 mins et un autre 4 mins/4 mins 30. La règle du pouce donne ds estimateurs très différents de la validation croisée et de la méthode de Sheather et Jones.