STAT2—Statistique non paramétrique École Centrale de Nantes

Feuille d’exercices

Exercice 1 (Propriétés statistique de Parzen–Rosenblatt).

On s’intéresse ici aux propriétés statistiques de l’estimateur de Parzen–Rosenblatt.
a) Rappelez l’expression de cet estimateur.
b) Montrez que si f est au moins C 2 alors pour tout x ∈ R

h2 00
Biais{fˆh (x)} = f (x)µ2 (K) + o(h2 ), h → 0,
2
R
où µ2 (K) = u2 K(u)du.
c) Montrez que pour K ∈ L2 , on a pour tout x ∈ R
Z  
ˆ 1 2 1
Var{fh (x)} = f (x) K(u) du + o , nh → ∞.
nh nh

d) En déduire une expression approchée pour la MISE.

e) Trouvez la fenêtre optimale minimisant cette MISE approchée.

Solution 1. a)
n  
1 X Xi − x
fˆh (x) = K , K noyau.
nh i=1 h

b) Puisque les Xi sont iid on a facilement

  
ˆ −1 X −x
E{fh (x)} = h E K
h
 
u−x
Z
−1
=h K f (u)du
h
Z
= K(ũ)f (x + hũ)dũ, ũ = (u − x)/h.

Faisons un développement de Taylor à l’ordre 2 pour f , on a

1
f (x + hũ) = f (x) + hũf 0 (x) + (hũ)2 f 00 (x) + o(h2 ).
2
De sorte que pour h → 0, on a
Z Z Z
ˆ 0 1 2 00
E{fh (x)} = f (x) K(u)du + hf (x) uK(u)du + h f (x) u2 K(u)du + o(h2 )
2
2
h
= f (x) + 0 + f 00 (x)µ2 (K) + o(h2 ).
2
c) Les Xi étant iid on a
  
ˆ −1 −2 X −x
V ar{fh (x)} = n h V ar K .
h

1
 De plus comme V ar(Y ) = E(Y 2 ) − E(Y )2 , il nous reste à calculer le premier terme (le 2ème
étant déjà fait). C’est parti pour le deuxième !
(  2 )  2
− u−x
Z
−1 −2 X x −1 −2
n h E K =n h K f (u)du
h h
Z
−1 −1
=n h K(u)2 f (x + hu)du
Z 
−1 −1 2
=n h K(u) f (x)du + o(1)
Z  
1 2 1
= f (x) K(u) du + o .
nh nh

Quand au premier on avait trouvé

  2
−1 −2 X −x
n h E K = n−1 h−2 {f (x) + o(h)}
h

qui est négligeable devant o(1/nh).

On a donc le résultat attendu.
d) On a donc puisque M SE = Biais2 + V ariance
Z 4 Z Z  
h 00 2 2 1 2 4 1
M ISE = {f (x)} µ2 (K) dx + f (x) K(u) dudx + o(h ) + o
4 nh nh
4 Z Z  
h 1 1
= µ2 (K)2 {f 00 (x)}2 dx + K(u)2 du + o(h4 ) + o .
4 nh nh

e) Reste plus qu’à dériver par rapport à h et résoudre (on néglige les petits o).
Z Z
3 2 00 1
h µ2 (K) {f (x)} dx − 2 K(u)2 du = 0
2
nh
Z Z
5 2 00 2 1
⇐⇒h µ2 (K) {f (x)} dx − K(u)2 du = 0
n
R
K(u)2 du
⇐⇒h5 = R
nµ2 (K)2 {f 00 (x)}2 dx
R 1/5
K(u)2 du


PPP
⇐⇒h = R .
nµ2 (K)2 {f 00 (x)}2 dx

Exercice 2 (La règle empirique de Silverman).

Dans cet exercice nous allons essayer de mieux comprendre ce qu’il se cache derrière la règle
empirique de Silverman et voir quelques modifications de cette dernière.
a) Soit f ∼ N (µ, σ 2 ). Montrez que
Z
3
f 00 (x)2 dx = √ .
8σ 5 π

Astuce : on utilisera le fait que pour Y ∼ N (µ, σ 2 ), E{(X − µ)2k } = (2k)!σ 2k /(2k k!).

2
b) En déduire que le choix par défaut de la fenêtre selon Silverman est donné par
v
 5 1/5 u n
4σ̂ u 1 X
hSilverman = , σ̂ = t (Xi − X̄)2
3n n − 1 i=1

c) Expliquez l’intuition derrière la formule suivante

 1/5
X[3n/4] − X[n/4]
 
4
h∗ = min σ̂, ,
3n 1.349
où X[np] représente la [np]–ième statistique d’ordre, i.e., la [np]–ième plus petite valeur de
l’échantillon X1 , . . . , Xn .
Solution 2. a) On commence par calculer f 00 (x). Rien de bien compliqué et on trouve
(x − µ)2
 
00 1
f (x) = − 2 f (x).
σ2 σ
Ensuite on notera que
1
f (x)2 = √ f (x; µ, σ 2 /2),
2 πσ 2

de sorte que
2
(x − µ)2
Z Z 
00 2 1
f (x) dx = − 2 f (x)2 dx
σ4 σ
(X − µ)4 (X − µ)2 σ2
      
1 1 2
= √ E + 4 − 2E , X∼N 0,
2 πσ 2 σ8 σ σ σ4 2
4!(σ/2)4 σ 2 /2
 
1 1
= √ + 4 −2 6
2 πσ 2 22 2!σ 8 σ σ
3
= 5√ .
8σ π
Nous avons utilisé pour l’avant dernière égalité le fait que
4! 4
E{(Y − µ)4 } = σ , Y ∼ N (µ, σ 2 ).
22 2!
b) D’après nos calculs précédents on sait déjà que (ϕ ∼ N (0, 1))
Z
1
ϕ(x)2 dx = √ .
2 π
De plus comme Z
µ2 (ϕ) := x2 ϕ(x)dx = 1,

on trouve alors que

Z 1/5 Z −1/5
−2/5 −1/5 2 002
h∗ = µ2 (ϕ) n ϕ f
−1/5
√

−1/5 3
=n (2 π)−1/5 4
√
8σ π
1/5
4σ 5

= .
3n
Comme σ est inconnue on l’estimera par sa version empirique σ̂.

3
c) Il suffit d’utiliser un estimateur robuste de la variance pour le cas Gaussien. Soit X ∼ N (µ, σ 2 ).
Alors puisque pour une loi N (0, 1) l’intervale inter-quartile vérifie

z0.75 − z0.25 = Φ(0.75) − Φ(0.25) ≈ 1.349,

et puisque (X − µ)/σ ∼ N (0, 1) ceci suggère d’estimer σ à l’aide de l’équation

X[3n/4] − µ X[n/4] − µ
− ≈ 1.349
σ σ
X[3n/4] − X[n/4]
⇐⇒ ≈ 1.349,
σ
soit l’estimateur suivant
X[3n/4] − X[n/4]
σ̃ = .
1.349
d) Vous savez que l’estimateur de l’écart-type σ
v
u n n
u 1 X 1X
σ̂ = t (Xi − X̄)2 , X̄ = Xi ,
n − 1 i=1 n i=1

est très sensible aux valeurs extrêmes présentes dans l’échantillon. Comment feriez vous pour
palier à ce problème ?
e) On remarque que h∗ correspond à la règle empirique de Silverman pour laquelle on estime
l’écart-type comme le minimum entre l’estimateur classique et sa version robuste. C’est une
manière de combiner les deux approches précédentes en prenant la fenêtre la plus petite des

PPP
deux afin de ne pas obtenir une densité trop lissée.

Exercice 3 (Validation croisée  Leave one out ).

On suppose que la vraie densité vérifie f ∈ L2 et on pose h > 0.
a) Rappelez l’expression de ce type de validation croisée pour l’estimateur de Parzen–Rosenblatt.
b) Montrez que Z
E{CV (h)} = M ISE(h) − f (x)2 dx.

c) Qu’en déduisez vous ?

Solution 3. a) L’expression est donnée par

n n  
Xj − Xi
Z
2 XX
CV (h) = fˆh (x)2 dx − K .
n(n − 1)h i=1 j=1 h
j6=i

4
b) Faisons les calculs séparément. D’une part on a
Z  Z n o
E fˆh (x) dx = E fˆh (x)2 dx
2

Z n o2 
= E fˆh (x) − f (x) + f (x) dx
Z Z Z n o
= M SE{fˆh (x)}dx + f (x) dx − 2 E fˆh (x) − f (x) f (x)dx
2

Z Z n o Z
ˆ ˆ
= M ISE(fh ) + f (x) dx + 2 E fh (x) f (x)dx − 2 f (x)2 dx
2

Z h n oi
ˆ 2 ˆ
= M ISE(fh ) − f (x) dx + 2EX E fh (X)
   
−
Z
2 X 1 X
= M ISE(fˆh ) − f (x) dx + EX E K
2
h h

D’autre part on a
 
 
 n X n     
2 X X j − Xi 2 X2 − X1
E K = E K .
n(n − 1)h  h  h h
 i=1 j=1
 

j6=i

En comibinant les deux résultats on trouve bien le résultat espéré.

c) En moyenne minimiser CV (h) revient à minimiser M ISE(fˆh ) puisque le terme f (x)2 dx ne
R

PPP
dépend pas de h.

Exercice 4 (Old faithful geyser).

Dans cet exercice nous allons mettre tout ce que nous avons vu sur l’estimation non paramétrique
d’une densité de probabilité en s’appuyant sur le jeu de données old faithfull geyser. Ce jeu
de données collecte (entre autre) le temps d’attente entre deux éruptions du geyser Old Faithful
situté dans le parc de Yellowstone.
a) Importez le jeu de données et renseignez vous sur ce dernier via les commandes R
data(faithful)
?faithful

data(faithful)

b) Lisez la documentation de la fonction density.

c) Exécutez les commandes suivantes, dites ce qu’elles font et commentez les résultats
par(mfrow = c(1, 3), mar = c(4, 5, 0.5, 0))
for (bandwidth in c(0.5, 10, 4)){
plot(density(faithful$waiting, kernel = "gaussian", bw = bandwidth),
main = "")
rug(faithful$waiting)
}

5
 PPP
Exercice 5 (Mélange de gaussiennes).
Soit la fonction  
0.7 x−1
f (x) = 0.3ϕ(x) + ϕ , x ∈ R,
0.3 0.3
où ϕ(·) correspond à la densité d’une N (0, 1).
a) Montrez que f est une densité de probabilité.
b) Ecrivez une fonction R qui génère un n–échantillon (iid) selon cette loi.
c) Simulez un n–échantillon (n choisi par vos soins) et obtenez une estimation de la densité. Vosu
choisirez une fenêtre  optimale à l’oeil .
d) Sur un même graphqiue, comparer cette estimation à la densité théorique.
R
Solution 4. a) Clairement f est positive. Reste à montrer que f (x)dx = 1. On a
Z Z
0.7
f (x)dx = 0.3 + ϕ(u)0.3du, u = (x − 1)/0.3
0.3
= 0.3 + 0.7 = 1,
R
où nous avons utilisé que ϕ = 1.

PPP
b) C’est juste du R.

Exercice 6 (Nadaraya–Watson).
Dans cet exercice, nous allons retrouver la forme de l’estimateur de Nadaraya–Watson pour la
régression non paramétrique.
a) Soit K1 et K2 deux noyaux sur R montrez que le noyau (x, y) 7→ K1 (x)K2 (y) est un noyau sur
R2 .
b) Considérons l’estimateur de la densité bivariée f (x, y) suivant
n    
1 X Xi − x Yi − y
fˆh1 ,h2 (x, y) = Kh1 Kh2 .
nh1 h2 i=1 h1 h2

Montrez que
n  
Xi − x
Z
1X
y fˆh1 ,h2 (x, y)dy = Kh1 Yi .
n i=1 h1

c) En déduire l’expression de l’estimateur de Nadaraya–Watson pour la régression non paramétrique.

d) Lisez la documentation de la fonction ksmooth et analysez le code suivant
data(faithful)
attach(faithful)
plot(eruptions, waiting)
fit <- ksmooth(eruptions, waiting, kernel = "normal")
lines(fit, col = "seagreen3", lwd = 2)

e) Jouez un peu avec l’argument bandwidth pour faire le lien avec le cours.
f) Ecrivez un bout de code R permettant de choisir une fenêtre adaptée par leave-one-out.

6
Solution 5. a) Il s’agit de vérifier la positivité, l’intégrabilité à 1 et la symétrie. C’est trivial.
b) On a
n    
− −
Z Z
1 X X i x Y i y
y fˆh1 ,h2 (x, y)dy = y Kh1 Kh2 dy
nh1 h2 i=1 h1 h2
n  Z
1 X Xi − x
= Kh1 (Yi − h2 ỹ)Kh2 (ỹ)(−h2 dỹ), ỹ = (Yi − y)/h2
nh1 h2 i=1 h1
n  
1 X Xi − x
= Kh1 Yi , symétrie du noyau et intégrale à 1.
nh1 i=1 h1

c) On utilise simplement la relation f (y | x) = f (x, y)/f (x). On estime alors le numérateur via
la question précédente et le dénominateur par l’estimateur de Parzen–Rosenblatt. On trouve
alors
n−1 ni=1 Kh1 {(Xi − x)/h1 }Yi
P
r̂(x) = −1 Pn ,
n i=1 Kh1 {(Xi − x)/h1 }
qui est bien l’expression donnée en cours.
d) Bon là c’est à vous de faire le travail.
e) L’argument bandwidth correspond à la fenêtre h du cours.
f) Voici mon implémentation du leave-one-out
bandwidths <- seq(0.1, 2, length = 100)
mse <- rep(NA, length(bandwidths))
n.obs <- nrow(faithful)
attach(faithful)

for (i in 1:length(bandwidths)){
h <- bandwidths[i]

residuals <- rep(NA, n.obs)

for (j in 1:n.obs){
pred <- ksmooth(eruptions[-j], waiting[-j], kernel = "normal", bandwidth = h,
x.points = eruptions[j])

residuals[j] <- pred$y - waiting[j]

}

mse[i] <- mean(residuals^2)

}

plot(bandwidths, mse, type = "l")

7
 36
35
mse

34
33

0.5 1.0 1.5 2.0

bandwidths

## Fenetre optimale
hopt <- bandwidths[which.min(mse)]

plot(eruptions, waiting)
lines(ksmooth(eruptions, waiting, kernel = "normal", bandwidth = hopt),
col = "seagreen3", lwd = 2)

8
 ●
●
● ●
●
●
90

● ●● ●●●
● ● ●
● ● ● ● ●●
● ●
● ●● ● ●
● ● ● ●●●
●● ●● ●●● ●● ●
● ● ● ● ● ●●● ● ●
● ● ●●●●
● ●● ● ● ●
● ● ● ●● ● ● ● ● ●● ●
80

● ● ● ● ●● ●
● ● ● ● ●● ● ● ●●
● ●
● ● ● ●● ● ● ●● ●● ●●
●● ●● ● ● ● ● ● ●
●● ●● ● ● ● ● ●
● ● ●
● ● ●
●●
● ● ● ● ●●
● ● ● ●● ●
●
waiting

● ● ● ● ●
70

● ● ●●
● ●
●
●
● ●
● ● ●
● ● ● ●
● ● ●
● ● ● ●
60

● ● ● ●● ●
● ●● ● ● ●●
● ●● ●
● ●
●● ●
● ●● ● ● ●
●
●●●●● ● ●
●● ●● ● ●
● ● ● ● ●
●●●● ● ●
50

● ● ● ●
●●● ● ●
● ● ●
● ● ●
●●
● ●
●● ●
●

1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0

eruptions

PPP