Stat Nonp P20 Annotations
Stat Nonp P20 Annotations
1 Introduction et rappels 4
1.1 Qu’est-ce que la statistique non-paramétrique ? . . . . . . . . . . . . . 4
1.2 Quelques problèmes de statistique non-paramétrique . . . . . . . . . . 5
1.2.1 Estimation de la fonction de répartition . . . . . . . . . . . . . 5
1.2.2 Estimation de densité . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.3 Régression non-paramétrique . . . . . . . . . . . . . . . . . . 5
1.2.4 Tests non-paramétriques . . . . . . . . . . . . . . . . . . . . . 6
1.2.5 Classification supervisée . . . . . . . . . . . . . . . . . . . . . 6
1.2.6 Classification non-supervisée, exemple génération . . . . . . . 7
1.3 Rappels d’inégalités classiques . . . . . . . . . . . . . . . . . . . . . . 7
1.3.1 Inégalité de Markov . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.2 Inégalité de Bienaymé-Tchebycheff (B-T) . . . . . . . . . . . . 7
1.3.3 Inégalité de Hoeffding . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Théorèmes de convergence classique . . . . . . . . . . . . . . . . . . . 8
1.4.1 Lemme de Slutsky . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2 Delta-méthode . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Petits rappels sur l’espérance conditionnelle . . . . . . . . . . . . . . 9
1.5.1 Calcul d’espérance conditionnelle . . . . . . . . . . . . . . . . 9
1.5.2 Propriété du transfert conditionnel . . . . . . . . . . . . . . . 10
1.6 Rappels sur les quantiles et les lois symétriques . . . . . . . . . . . . 11
1.6.1 Quantiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6.2 Loi symétrique . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.7 Rappels sur les tests (cadre paramétrique) . . . . . . . . . . . . . . . 12
1.7.1 Comparaison de test, principe de Neyman . . . . . . . . . . . 15
1.7.2 Explications sur des exemples . . . . . . . . . . . . . . . . . . 16
1.7.3 La p-valeur . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.7.4 Interprétation des p-valeurs : d’autres exemples et détails . . . 26
1.8 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1
2.4 Test d’homogénéité de Kolmogorov Smirnov . . . . . . . . . . . . . . 43
2.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3 Tests robustes 52
3.1 Un exemple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2 Un test paramétrique : le test de Student . . . . . . . . . . . . . . . . 53
3.2.1 Un seul échantillon . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2.2 Deux échantillons indépendants . . . . . . . . . . . . . . . . . 54
3.2.3 Echantillons appariés (paired data) . . . . . . . . . . . . . . . 55
3.2.4 Importance des conditions d’application . . . . . . . . . . . . 57
3.3 Test du signe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.3.1 Test du signe sur un seul échantillon . . . . . . . . . . . . . . 61
3.3.2 Test du signe sur deux échantillons . . . . . . . . . . . . . . . 63
3.4 Statistiques d’ordre et de rang . . . . . . . . . . . . . . . . . . . . . . 66
3.5 Test des rangs signés de Wilcoxon . . . . . . . . . . . . . . . . . . . . 66
3.5.1 Sur un échantillon . . . . . . . . . . . . . . . . . . . . . . . . 66
3.5.2 Echantillons appariées . . . . . . . . . . . . . . . . . . . . . . 73
3.6 Wilcoxon de la somme des rangs/Mann-Whitney . . . . . . . . . . . 74
3.6.1 Résultats préliminaires sur le vecteur des rangs . . . . . . . . 74
3.6.2 Test de Mann-Whitney . . . . . . . . . . . . . . . . . . . . . . 76
2
Introduction
Ces notes de cours font suite aux notes du cours d’introduction à la statistique
non paramétrique de Catherine Mathias, Vincent Rivoirard et Laëtitia Comminges.
3
Chapitre 1
Introduction et rappels
Yi = f (Xi ) + i , i = 1, ..., n
4
1.2 Quelques problèmes de statistique non-paramétrique
1.2.1 Estimation de la fonction de répartition
On observe X1 , . . . , Xn n variables réelles de loi P . On cherche à estimer la loi
P . Or P est entièrement décrite par sa fonction de répartition
R → [0, 1]
F :
x → P (] − ∞, x])
Yi = f (Xi ) + i , i = 1, . . . , n
5
1.2.4 Tests non-paramétriques
Deux exemples de problèmes possibles :
— Soit X une v.a. et P une distribution donnée. A l’aide de X1 , . . . , Xn iid de
même loi que X, tester :
H0 : X ∼ P, contre X 6∼ P
6
1.2.6 Classification non-supervisée, exemple génération
Idée : ayant quelques réalisations X1 , ..., Xn i.i.d. de la loi P comment générer
d’autres instances de Y1 , Y2 suivant la même loi ? Exemple : peintures de paysages.
Ce sont des algorithmes de type GAN (Generative Adversarial Networks), VAE
(Variational auto-encoder), ...
E(X)
P(X ≥ t) ≤
t
Var(X)
P(|X − E(X)| > t) ≤
t2
7
1.3.3 Inégalité de Hoeffding
Soient Y1 , . . . , Yn des v.a.r. indépendantes centrées et telles que
Alors n
X 2t2
∀t > 0, P( Yi ≥ t) ≤ exp − Pn
i=1 (bi − ai )
2
i=1
iid
Remarque 1.4. Comparaison entre Hoeffding et B-T : soit X1 , . . . , Xn ∼ Be(p)
avec p ∈ (0, 1). On cherche un intervalle de confiance bilatéral à gauche de niveau
1 − α pour p avec l’une des inégalités ci-dessus :
— Si on utilise l’inégalité B-T, P(|X̄ − p| > c) ≤ p(1−p)
nc2
≤ 4c12 n := α. Donc
P(p ∈ [X̄ − 2√1nα , X̄ + 2√1nα ]) ≥ 1 − α. Pour α = 5% et n = 100, la précision,
i.e. la longueur, de cet intervalle est √1nα = 0.22.
— Si on utilise l’inégalité
r Hoeffding, P(|X̄ − p| > c) ≤ 2 exp(−2nc2 ), Donc
de r
2 2
log( α ) log( α )
P(p ∈ [X̄ − 2n
, X̄ + 2n
]) ≥ 1 − α. Pour α = 5% et n = 100, la
r
2
log( α )
précision, i.e. la longueur, de cet intervalle est 2 n
= 0.14.
1.4.2 Delta-méthode
On se donne une suite (Un )n de vecteurs aléatoires de Rm , une suite déterministe
(an )n et une application ` : Rm → Rp telles que
— an → +∞
— ∃U ∈ Rm un vecteur déterministe (=constant) et V un vecteur aléatoire tels
loi
que an (Un − U ) → V .
8
— ` est une fonction différentiable en U de différentielle D`(U ) ∈ Mpm (R).
Alors on a la convergence en loi
loi
an (`(Un ) − `(U )) → D`(U )V.
iid
Exemple 1.5. Soit X1 , . . . , Xn ∼ P (λ) avec λ > 0. Alors d’après le TCL on a
√ loi
n(X̄ − λ) → N (0, λ). Donc on a aussi, d’après le théorème ci-dessus,
√ q √ loi 1
n( X̄ − λ) → N (0, )
4
√ √
En effet ici Rm = Rp = R, Un = X̄, U = λ, an = n, V ∼ N (0, λ), `(u) = u,
donc D`(U ) = `0 (λ) = 2√1 λ et D`(U )V ∼ N (0, 4λ
λ
).
Exemples :
1. soient Z et T deux variables aléatoires indépendantes de loi exponentielle de
paramètre λ. On note S = Z + T et on cherche à calculer la variables aléatoire
E(Z | S). Soit s > 0. On trouve facilement (car tout le monde a une densité ...
) que la densité conditionnelle de Z sachant S = s est Rs
donnée par fZS=s (z) =
1
1 (z). On a alors immédiatement E(Z | S = s) = 0 zfZS=s (z)dz = 2s . Et en
s [0,s]
utilisant la propriété que l’on vient de rappeler, on a finalement E(Z | S) = S2 .
2. Soit U et V deux v.a. réelles.
h
On rappelle la définition
i
de la variance
condi-
2
tionnelle : Var(U | V ) = E (U − E(U | V ))2 | V = E[U 2 | V ] − E[U | V ] .
On a h i
E Var(U | V ) = E[g(V )] avec g(v) = Var(U | V = v).
2
En effet Var(U | V ) = E[U 2 | V ] − E[U | V ] = `(V ) − (h(V ))2 avec `(v) =
E[U 2 | V = v] et h(v) = E[U | V = v]. Donc `(v) − h(v)2 = Var(U | V = v).
9
1.5.2 Propriété du transfert conditionnel
Soit f : Rp+k → Rq une fonction borélienne. Alors la loi conditionnelle de f (X, Y )
sachant X = x vérifie
PX=x X=x
f (X,Y ) = Pf (x,Y )
et donc h i h i
E f (X, Y ) | X = x = E f (x, Y ) | X = x
En particulier si X et Y sont indépendantes, on a
PX=x
f (X,Y ) = Pf (x,Y )
et donc h i h i
E f (X, Y ) | X = x = E f (x, Y ) .
où
g(x) = E[1Y ≤X | X = x] = E[1Y ≤x ] = FY (x)
où on a noté FY la cdf de Y . Donc on a
iid
Exemple : Reprenons l’exemple de la sous-section précédente : Z, T ∼ exp(λ). On
veut calculer E(S 2 Z | S = s) pour s > 0. En utilisant la propriété ci-dessus on
3
obtient E(S 2 Z | S = s) = E(s2 Z | S = s) = s2 .
10
1.6 Rappels sur les quantiles et les lois symétriques
1.6.1 Quantiles
On ne donne ici que la définition dans le cas simple où la loi est de cdf F continue
et strictement croissante.
Soit X une variable aléatoire réelle de cdf F continue et strictement croissante.
Pour α ∈ (0, 1), on appelle quantile d’ordre α de la loi F l’unique réel qαF tel que
F (qαF ) = P (X ≤ qαF ) = α
autrement dit
qαF = F −1 (α) (1.1)
Attention, quand la cdf n’est pas continue, l’équation ci-dessus n’a pas toujours
de solution. De plus si la cdf n’est pas strictement croissante, l’équation peut avoir
une infinité de solutions. La définition générale d’un quantile sera vue dans le cha-
pitre 2.
et
P(X > 0) = P(X < 0) (1.3)
(1.2) se réécrit
11
ce qui implique
H0 : θ ∈ Θ0
contre « l’alternative »
H1 : θ ∈ Θ1 ,
avec Θ0 ⊂ Θ, Θ1 ⊂ Θ et Θ0 ∩ Θ1 = ∅. Construire un test signifie construire une
procédure φ = φ(X) de la forme
0 si X ∈
/ R. « on accepte l’hypothèse nulle »
φ(X) = 1{X∈R} = (1.6)
1 si X ∈ R. « on rejette l’hypothèse nulle »
avec R mesurable.
12
o 1.7. On désigne indifféremment l’ensemble R ⊂ A ou bien l’événement
Définition
n
X ∈ R comme zone de rejet ou encore zone critique du test φ.
Erreur de test
Lorsque l’on effectue un test, il y a quatre possibilités. Deux sont anecdotiques
et correspondent à une bonne décision :
— Accepter l’hypothèse H0 alors que θ ∈ Θ0 (c’est-à-dire l’hypothèse H0 est
vraie).
— Rejeter l’hypothèse H0 alors que θ ∈ Θ1 (c’est-à-dire l’hypothèse H0 est
fausse).
Les deux autres possibilités sont celles qui vont nous occuper, et correspondent
à une erreur de décision :
— Rejeter l’hypothèse H0 alors que θ ∈ Θ0 (c’est-à-dire l’hypothèse H0 est vraie).
— Accepter l’hypothèse H0 alors que θ ∈ Θ1 (c’est-à-dire l’hypothèse H0 est
fausse).
Définition 1.9. [Erreur de première et seconde espèce] L’erreur de pre-
mière espèce, ou encore de "type I" correspond à la probabilité maximale
de rejeter l’hypothèse alors qu’elle est vraie :
h i h i
sup Eθ φ(X) = sup Pθ X ∈ R .
θ∈Θ0 θ∈Θ0
13
de première espèce revient à faire un « faux négatif », et commettre une
erreur de seconde espèce revient à faire un « faux positif ».
β : Θ1 → [0, 1]
définie par h i
θ ∈ Θ1 ; β(θ) = Pθ X ∈ R .
Une illustration intuitive des erreurs et paramètres α et β est donnée en figure 1.4.
14
1.7.1 Comparaison de test, principe de Neyman
Idéalement, on souhaite que l’erreur de première espèce et l’erreur de seconde
espèce soient toutes deux simultanément petites. Les deux tests triviaux
φ1 = 1∅ , et φ2 = 1X
Remarque 1.12. On ne peut pas toujours faire en sorte que l’erreur de première
espèce soit égale à α (problème de non continuité d’une fonction de répartition par
exemple, cf chapitre 2 en particulier). C’est pourquoi on se contente d’exiger que
l’erreur de première espèce soit plus petite que α.
Définition 1.13. On dit qu’un test est de taille α si l’erreur de première espèce
est égale à α.
15
iid
On observe donc X1 , . . . , X20 ∼ N (µ, σ 2 ) avec µ et σ 2 inconnus. Pour le directeur
de l’usine, l’erreur la plus grave serait de conclure que le niveau de polluant est trop
élevé alors qu’il ne l’est pas. Il choisit donc comme hypothèses
H0 : µ ≤ 6 contre H1 : µ > 6.
Prenons maintenant le point de vue de l’écologiste. Si la limite est supérieure à
8mg/kg, il y a danger. Contrairement au directeur d’usine, l’écologiste considère
que l’erreur la plus grave serait de conclure que le niveau de polluant n’est pas trop
élevé alors qu’en réalité il l’est. Il effectue donc le test suivant
H0 : µ ≥ 8 contre µ < 8.
La mise en oeuvre de ces tests sera faite en exercice (cf TD1 exercice 2).
H0 : µ = 3 contre µ 6= 3
Dans la pratique : nous observons x1 , . . . , xn . Comme nous voulons savoir si la
moyenne est égale à 3 ou plus grande que 3, naturellement nous regardons la moyenne
empirique x̄. Imaginons que x̄ = 3.5. Alors que conclure ? Et bien ça dépend...ça
dépend de plusieurs facteurs, plus exactement, ça dépend ici de n et de σ.
En effet, le problème est que, évidemment, on ne tombera jamais sur 3 exacte-
ment. Imaginons que la vraie moyenne est 3. Alors comme les Xi sont aléatoires et
qu’on n’en a qu’une quantité finie n, on n’a jamais l’information exacte sur µ en
utilisant l’échantillon, mais seulement une information approchée et aléatoire.
Donc si la moyenne empirique vaut 3.5, la question est : est-ce que la vraie
moyenne est 3 et que je tombe sur 3.5 parce que c’est aléatoire ? Ou bien est-ce que
c’est parce que ce n’est pas 3 la vraie moyenne ?
Pour répondre à ces questions, il faut utiliser les tests, et surtout utiliser toutes
les informations que l’on a à notre disposition (ou que l’on peut déduire des données),
en particulier la taille de l’échantillon et la variance σ 2 . En effet ce sont ces deux
informations qui vont nous aider à savoir si c’est "normal" de tomber sur 3.5 en
ayant une vraie moyenne de 3, ou bien si c’est "anormal" (ou "atypique").
Ici regardons ce qui se passe sous H0 , c’est-à-dire quand µ = 3 (c’est toujours
ce qu’on fait en fréquentiste, on regarde ce qu’il est censé se passer sous H0 , donc
asymétrie des deux hypothèses). Si on est vraiment sous H0 , alors la question est :
qu’est-ce qu’une valeur "normale" (ou "usuelle" ou "typique") de X̄ quand µ = 3 ? Il
suffit pour cela de standardiser pour se ramener à une variable normale standard et
utiliser les quantiles de N (0, 1). En effet, si µ = 3 on a
√ X̄ − 3
n ∼ N (0, 1)
σ
Donc, comme le quantile d’ordre 97.5% de la loi normale standard vaut 1.96 (environ)
on a
√ |X̄ − 3|
P3 n ≤ 1.96 = 95%
σ
16
Autrement dit, on peut dire que, avec une très grande probabilité, ici plus précisé-
ment avec une probabilité de 95%, la variable aléatoire
√ X̄ − 3
T = n
σ
se trouve dans l’intervalle [−1.96, 1.96]. Autrement dit, une valeur "typique" de la
statistique T , si on est vraiment sous H0 , est une valeur entre -1.96 et 1.96.
Ainsi si on tombe sur une valeur qui sort de cette intervalle, on se dit que ça
n’est pas une valeur "normale" pour T sous H0 et donc on rejette H0 .
Il est évidemment toujours possible que, tout en étant sous H0 , c’est-à-dire ici,
tout en ayant une vraie valeur de µ égale à 3, on tombe sur une valeur observée de
T qui sorte de l’intervalle [−1.96, 1.96], puisque la loi normale a son support sur R.
Mais ceci se produit "rarement" et donc la possibilité de se tromper en rejetant à
tort H0 est faible : ici 5% (on prend toujours α petit). C’est l’erreur de type I.
Maintenant illustrons cette dépendance par rapport à σ et n dans notre exemple
(donc on suppose toujours x̄ = 3.5) sur notre décision finale .
1. Imaginons
√ d’abord que σ = 1 et n = 100. Alors la valeur observée de T est
3.5−3
t = 100 1 = 5. Comme 5 est en dehors de l’intervalle [−1.96, 1.96], on
conclut que c’est une valeur "anormale" pour H0 et donc on rejette H0 .
2. √Imaginons que σ = 5 et n = 100. Alors la valeur observée de T est t =
100 3.5−3
5
= 1. Alors on accepte H0 . L’idée est que c’est très possible que la
vraie valeur de µ soit 3 et de tomber sur une valeur aussi grande que 3.5 ici,
car les données ont une grande variance.
√
3. Imaginons que σ = 1 et n = 9. Alors la valeur observée de T est t = 9 3.5−3 1
=
1.5. Alors on accepte à nouveau H0 . L’idée est qu’une valeur de x̄ = 3.5 n’est
pas "anormale" pour H0 si on n’a pas beaucoup de données (le résultat est
peu précis si on n’a très peu de données donc il n’est pas "anormal" d’avoir
vraiment µ = 3 tout en ayant une valeur x̄ un peu "éloignée" de 3).
17
C’est-à-dire que, avec une grande probabilité, plus précisément ici 95%, et si on est
vraiment sous H0 , la statistique T doit être plus petite que 1.64. Donc on rejette H0
si ce n’est pas le cas.
Fonction puissance
Maintenant, après avoir construit le test, on est intéressé par l’erreur de seconde
espèce et par la fonction puissance, c’est-à-dire, on est intéressé par ce qui se passe
sous H1 . On veut que la probabilité de rejeter H0 , quand on est sous H1 , soit grande,
c’est-à-dire qu’on veut que la puissance soit grande. Éventuellement la fonction
puissance nous permet de comparer différents tests. Une des propriétés souhaitées
est alors que, si on a suffisamment de données, on puisse dire qu’on est sous H1
quand on l’est bien, avec une très grande probabilité. C’est le cas quand le test est
"convergent" (ou "consistant") : la fonction puissance tend vers 1 quand n tend vers
l’infini.
Évidemment, comme son nom l’indique, la puissance est une fonction, car elle
dépend de l’alternative exacte. En effet en général, Θ1 est une hypothèse composite,
c’est-à-dire que Θ1 n’est pas un singleton et on a souvent une infinité de cas possibles
( exemple : Θ1 = R \ {3} ou Θ1 =]3, +∞[)). Il est évidemment plus facile de voir
qu’on est sous H1 quand le vrai µ vaut 10 que quand il vaut 3.5 (toutes choses
étant égales par ailleurs). De plus, la puissance dépend également de la taille de
l’échantillon et de sigma.
18
est donnée par
√ X̄ − 3
β(µ) = Pµ (|T | > q) = Pµ (| n | > q)
σ
√ X̄ − 3 √ X̄ − 3
= Pµ ( n > q) + Pµ ( n < −q)
σ σ
√ X̄ − µ + µ − 3 √ X̄ − µ + µ − 3
= Pµ ( n > q) + Pµ ( n < −q)
σ σ
√ X̄ − µ √ 3−µ √ X̄ − µ √ 3−µ
= Pµ ( n >q+ n ) + Pµ ( n < −q + n )
σ σ σ σ
√ X̄ − µ √ 3−µ √ X̄ − µ √ 3−µ
= 1 − Pµ ( n ≤q+ n ) + Pµ ( n < −q + n )
σ σ σ σ
√ 3−µ √ 3−µ
= 1 − Φ(q + n ) + Φ(−q + n )
σ σ
Quelques exemples d’applications numériques avec α = 5% (arrondis à deux
chiffres après la virgule) :
Code python pour calculer α
def calcul_puissance(alpha,sigma,n,mureel,muH0):
q = stat.norm.ppf(1.0-alpha/2,loc=0,scale=1)
beta=(1.0- stat.norm.cdf(q + np.sqrt(n)*(muH0-mureel)/sigma,loc=0,scale=1)
+ stat.norm.cdf(- q + np.sqrt(n)*(muH0-mureel)/sigma,loc=0,scale=1))
print("mu=",mureel," muH0=",muH0," sigma=",sigma," n=",n,
" beta(",mureel,")=",np.round(beta,2))
return beta
calcul_puissance(0.05,1,100,3.5,3.0);
calcul_puissance(0.05,1,10,3.5,3.0);
calcul_puissance(0.05,2,100,3.5,3.0);
calcul_puissance(0.05,1,100,3.1,3.0);
murange = np.linspace(0,6,100)
betan100=np.zeros_like(murange)
betan10=np.zeros_like(murange)
betan2=np.zeros_like(murange)
plt.figure(14)
plt.rc(’font’,size=14)
19
plt.plot(murange,betan100,"g",murange,betan10,"b",
murange,betan2,"r",linewidth=4)
plt.ylabel("Puissance",size=14)
plt.xlabel("$\mu$",size=14)
plt.legend(["n=100","n=10","n=2"])
plt.title("Puissance pour $\sigma=1, \mu_{H0}=3.0$")
plt.savefig("betaplot.pdf")
==============Resultats:==============================
0.8
Puissance
0.6
0.4
0.2
n=100
n=10
n=2
0 1 2 3 4 5 6
Figure 1.5 – Fonction puissance pour les exemples du test bilatéral µ = 3 contre
µ 6= 3.
√ X̄ − 3
β(µ) = Pµ (T > q) = Pµ ( n > q)
σ
√ X̄ − µ √ 3−µ
= Pµ ( n >q+ n )
σ σ
√ X̄ − µ √ 3−µ
= 1 − Pµ ( n ≤q+ n )
σ σ
√ 3−µ
= 1 − Φ(q + n )
σ
20
Dans ces deux premiers exemples, on voit immédiatement que la fonction puis-
sance tend vers 1 lorsque n → ∞. Cela signifie que pour tout µ de l’alternative
et pour tout > 0, il existe une taille d’échantillon n0 telle que la probabilité
de rejeter à tort H1 , quand on est sous Pµ pour ce µ particulier, est plus pe-
tite que si n ≥ n0 . En revanche, dans les deux cas, on peut montrer que la
fonction puissance ne tend pas vers 1 uniformément, ce qui signifie que ce n0
dépend de µ (considérer par exemple la suite µn = 3 + n1 ). L’erreur de seconde
espèce, qui est définie par un sup, ne tend pas vers 0. Voir aussi l’encadré 1.7.1.
3. σ inconnu et problème de test H0 : µ = 3 contre H1 : µ > 3.
√
Si σ est inconnu, on ne peut plus baser notre test sur T = n X̄−3 σ
car T
n’est
q plus calculable. On remplace donc σ par un estimateur, ici prenons σ̂ =
1 Pn √ X̄−3
i=1 (Xi − X̄) . Avec cet estimateur σ̂ on définit donc T =
2 n σ̂ . La
n−1
loi de cette statistique sous H0 est la loi de Student à n − 1 degrés de liberté.
En effet on a
√ X̄−3 √ X̄−3
√ X̄ − 3 n σ n σ
T = n = σ̂ = q (1.8)
σ̂ σ
σ̂ 2
2σ
√ √
qσ̂+ n(3−µ)
Si on pose U = n X̄−µσ
et V = σ
on a β(µ) = 1−P(U ≤ V ) avec U et
V indépendantes, puisque X̄ et σ̂ sont indépendantes. Donc, d’après√l’exemple
2 de la section 1.5, β(µ) = 1 − E[FU (V )] où FU est la cdf de U = n X̄−µσ
∼
N (0, 1). On obtient donc finalement
√
qσ̂ + n(3 − µ)
β(µ) = 1 − E Φ( ) .
σ
On peut ensuite vérifier si la puissance tend bien simplement vers 1 quand
n tend vers l’infini. Pour cela on peut utiliser l’expression ci-dessus. Mais
on peut aussi déduire cette propriété directement, sans faire appel à cette
21
expression. Rappelons une méthode assez fréquemment utilisée pour montrer
cette propriété : pour fixer les idées, on doit montrer que Pµ (T > cn ) tend
vers 1 quand n tend vers l’infini. Il suffit alors de :
— montrer que la statistique de test Tn se décompose en Tn = Tn,0 + Tn,1
avec
— Tn,0 = OP (1) ("grand O en probabilité", typiquement on montre que Tn,0
converge en loi)
P roba
— Tn,1 → +∞,
— et cn = O(1).
√ √ √
On a ici : Tn = n X̄−3
σ̂
= Tn,0 + Tn,1 avec Tn,0 = n X̄−µ
σ̂
, Tn,1 = n µ−3
σ̂
et
T (n−1)
cn = q1−α . On a, sous Pµ avec µ > 3 et quand n → +∞,
Tn,0 ∼ T (n − 1) donc Tn,0 = OP (1).
p.s.
Tn,1 → +∞
T (n−1) N (0,1)
et enfin q1−α → q1−α car la loi de Student à n − 1 degrés de liberté tend
vers la loi normale standard (cf chapitre 2 théorème 2.20).
1.7.3 La p-valeur
Définition 1.14. Supposons avoir construit une famille de tests φα (X), chacun de
niveau α, pour α ∈ [0, 1]. La p-valeur associée à cette famille est la variable aléatoire
réelle définie par
p(X) = inf{α ∈ [0, 1] : φα (X) = 1}.
Remarque 1.15. On constate que p(X) est le niveau à partir duquel on se met
à rejeter H0 . C’est comme si on faisait le test sans connaître le α et on tire la
conclusion à la fin une fois que le α est dévoilé. Donc
— Si p(x) < α alors on rejette H0 au niveau α.
— Si p(x) > α alors on conserve H0 au niveau α.
Mise en garde 1.7.2. Une p-valeur petite ne veut pas dire que
l’on a plus de chances d’être sous H1 que sous H0 : ça dépend en
fait du comportement de la p-valeur sous H1 . On sait que, sous cer-
taines conditions du moins (cf chapitre 2), la p-valeur suit une loi
uniforme sous H0 , mais on ne sait pas forcément le comportement
22
de la p-valeur sous H1 . La question de la probabilité de H0 sachant
les données est une question bayésienne à laquelle on peut répondre
si on a un a priori sur l’alternative (il faut aussi parfois un a priori
sur H0 ). Attention donc à l’interprétation des p-valeurs, ne pas
dire "la p-valeur est petite donc la probabilité que H0 soit fausse
est grande".
Pour autant, une p-valeur importante n’implique pas forcément que
H0 soit vraie. Il se peut que le test ne soit pas puissant. Par exemple
considérons le test φ(X) ≡ 0 : ce test accepte toujours H0 . L’en-
semble dans la définition de la p-valeur est vide, par convention on
prend son sup pour définir la p-valeur, c’est-à-dire que la p-valeur
est égale à 1.
p(x) = 1 − F0 (T (x)).
En effet,
α1 ≤ α2 =⇒ Rα1 ⊂ Rα2
23
autrement dit, si on rejette à un niveau α1 alors on rejette aussi à tout
niveau α2 ≥ α1 .
Admis.
Ce qu’on veut dire par loi fixe : T (X) a la même loi ∀θ ∈ Θ0 . Par exemple c’est
le cas si Θ0 = {θ0 }. C’est aussi le cas pour le test de Kolmogorov-Smirnov, le test
du signe et les tests de Wilcoxon (cf chapitre 2).
H0 : µ = 2 contre µ ≤ 2.
Le test utilisé est alors le test de Student (cf "exemples de calcul de puissance",
√
item 3). Le test est alors φ = 1T ≤qT (n−1) avec T = n X̄−2 σ̂
. On est alors dans
α
les conditions d’application du théorème de Wassermann, item 1 : en effet,
24
puisque la loi de Student est une loi continue, on a bien un test de taille α (et
pas seulement de niveau α). La p-valeur observée est donc donnée par
p(xn1 ) = FT (n−1) (T (xn1 ))
où T (xn1 ) est l’observation de la statistique T (X1n ).
P q P
1 20 1 20
i=1 (xi − x̄) =
Application numérique : n = 20, 20 2
i=1 xi = 1.34 et
√ 20
1.06, ce qui donne la valeur observée T (xn1 ) = 20(1.34 − 2)/1.06 = −2.78. La
p-valeur observée est donnée par p(xn1 ) = 0.01. On peut trouver cette valeur
sur R en utilisant la commande pt(-2.78,19). Pour obtenir directement ce
résultat sans faire de calcul, on peut utiliser la commande R
t.test(-2.78,mu=2,alternative="less").
iid
— Soit X1 , . . . , Xn ∼ Be(p). On veut tester
H0 : p = 1/2 contre p > 1/2.
P
On utilise la statistique T = ni=1 Xi qui suit, sous H0 , une loi binomiale
B(n, 1/2). Au vu de H1 , on rejette quand T est trop grand. Donc on pose
φ = 1T >c où c est déterminé par le fait que le test est de niveau α
P1/2 (T > c) ≤ α. (1.9)
Attention ici la statistique T est discrète donc sa cdf FB(n,1/2) sous P1/2 n’est
pas continue. Donc on ne peut pas toujours avoir l’égalité.
On verra au chapitre 2 que le plus petit entier c vérifiant (1.9) est donné par
B(n,1/2)
c = q1−α , i.e. le quantile d’ordre 1 − α de la loi binomiale B(n, 1/2). Cela
donne le test suivant
φ = 1T >qB(n,1/2) .
1−α
25
Méthode pour construire un test
1. Choix de H0 et H1 .
2. Détermination de T (X), la statistique de test. On doit connaitre sa loi sous
H0 . Evidemment on souhaite aussi que cette statistique ait un comportement
différent sous H0 et sous H1 pour pouvoir discriminer les deux hypothèses.
3. Allure de la zone de rejet en fonction de H1 (i.e. en fonction du comportement
de T (X) sous H1 ).
4. Observation de la réalisation T (x) de T (X).
5. Calcul de la p-valeur associée p(x) et comparaison à un seuil fixé par un non-
statisticien.
6. Conservation ou non de H0 .
26
data: x
t = -1.9561, df = 19, p-value =
0.06532
alternative hypothesis: true mean is not equal to 1.5
95 percent confidence interval:
0.6181763 1.5298237
sample estimates:
mean of x
1.074
On tombe sur une p-valeur d’environ 0.06 donc on accepte (tout juste) H0 au
niveau 5%. Là encore, comme la p-valeur est proche du niveau, on n’a pas une
confiance énorme en le résultat final.
Interprétation à l’aide du théorème de Wassermann
Reprenons un des exemples précédents : premier item de "exemples de calculs
de p-valeurs". La p-valeur s’écrit dans cet exemple p(xn1 ) = FT (n−1) (T (xn1 )). Dans
l’application numérique, la valeur observée de la statistique T est t = −2.78. Si on
est vraiment sous H0 , t est alors censée être la valeur observée d’une statistique
qui suit une loi de Student à 19 degrés de liberté, et la p-valeur mesure alors la
probabilité qu’une variable de Student à 19 degrés de liberté soit plus petite que
-2.78, c’est-à-dire la probabilité d’observer une valeur de T plus petite que -2.78
si on est vraiment sous H0 . Donc la p-valeur mesure en quelque sorte le côté
atypique de la valeur observée, par rapport à ce qu’il est censé se passer sous H0 .
Ici, si on était vraiment sous H0 , il y aurait une probabilité de 1% d’observer
une valeur inférieure ou égale à −2.78 pour la statistique T . Donc -2.78 est plutôt
une valeur atypique pour H0 et on penche donc pour le rejet de H0 .
27
1.8 Exercices
Exercice 1.1. Soit X une variable aléatoire réelle, absolument continue de densité
continue f , de fonction de répartition F . On observe un n-échantillon iid (X1 , . . . , Xn )
de même loi que X. On considère la statistique T qui ordonne l’échantillon dans le
sens croissant :
T (X1 , . . . , Xn ) = (X(1) , . . . , X(n) ),
avec X(1) ≤ X(2) ≤ · · · ≤ X(n) . (X(1) , . . . , X(n) ) s’appelle la statistique d’ordre.
1. On suppose pour cette question uniquement que les Xi sont seulement indé-
pendants et de lois continues (c’est-à-dire que les Xi sont indépendants et ont
tous une fonction de répartition Fi continue, mais pas forcément absolument
continue). Montrer que
P (∃ i 6= j : Xi = Xj ) = 0,
(a) Montrer que nYn converge en loi vers une loi exponentielle.
(b) Étudier la convergence en loi de Zn , puis sa convergence en probabilité et
L1 .
(c) Soit > 0. Calculer P[|Zn − 1| > ]. En déduire que Zn converge presque
sûrement.
(d) Rappeler les implications logiques entre les modes de convergence étudiés :
en loi, en probabilité, en norme L1 , en norme L2 , presque sûre.
H0 : µ ≤ 6 contre H1 : µ > 6.
28
1 P
20
2. On calcule à partir de ces données x̄ = 7 et σ̂ 2 = 19 2 2
i=1 (xi − x̄) = 2.4 .
Calculer la p-valeur observée et conclure si on choisit le niveau α = 5%.
iid
Exercice 1.3. On dispose d’un échantillon de loi Bernoulli de paramètre p : X1 , . . . , Xn ∼
Be(p).
1. Proposer une procédure de test pour le problème suivant
29
Chapitre 2
Estimation de la fonction de
répartition
Remarque 2.2. Pour insister sur le caractère aléatoire de F̂n , on peut écrire parfois
F̂n (ω, x) au lieu de F̂n (x). F̂n (ω, x) désigne donc la valeur de la cdf F̂n en x quand
l’observation est ω.
Remarque 2.3. On construit facilement F̂n car c’est une fonction en escalier.
P
Fixons ω et écrivons (x1 , . . . , xn ) = (X1 (ω), . . . , Xn (ω)). Alors F̂n (ω, x) = n1 ni=1 1xi ≤x
est la fonction de répartition de la variables aléatoire Z à valeurs dans {x1 , . . . , xn }
et telle que P(Z = xi ) = nk si la valeur xi apparait k fois dans {x1 , . . . , xn }. Par
exemple, si tous les xi distincts, alors F̂n (ω, ·) est la cdf de la loi uniforme sur
{x1 , . . . , xn }.
Soit (X(1) , . . . , X(n) ) la statistique d’ordre associée à X1n . On rappelle que cela
signifie que {X(1) , . . . , X(n) } = {X1 , . . . , Xn } et
30
La fonction F̂n (ω, ·) est discontinue aux points X(j) (ω). Elle a un saut égal au nombre
de fois où la valeur Xi (ω) apparait dans {X1 (ω), . . . , Xn (ω)}. En particulier si tous
les Xi (ω) sont distincts, i.e. X(j) (ω) < X(j+1) (ω) pour tout j, alors F̂n (ω, x) = nj
pour tout x ∈ [X(j) (ω), X(j+1) (ω)[. Dans tous les cas, elle vaut 0 sur ] − ∞, X(1) (ω)[
et 1 sur [X(n) (ω), +∞[.
Proposition 2.4. Soit x ∈ R, F̂n (x) est un estimateur sans biais de F (x)
et limn→∞ F̂n (x) = F (x) p.s. Par ailleurs
√ Loi
n(F̂n (x) − F (x)) −→ N (0, F (x)(1 − F (x)))
P iid
Démonstration. F̂n (x) = n1 ni=1 1Xi ≤x avec 1Xi ≤x ∼ Be(F (x)) donc limn→∞ F̂n (x) =
F (x) p.s. découle de la LGN. La deuxième propriété vient du théorème limite central
en remarquant que Var(1Xi ≤x ) = F (x)(1 − F (x)).
Ce résultat est de nature paramétrique car x est fixé. On peut aller plus loin.
Démonstration. 1. Par définition de F (−1) (q), il existe une suite (un )n≥0 telle que
F (un ) ≥ q et un → F (−1) (q) en décroissant (c’est donc une limite à droite) .
n→∞
Comme F est continue à droite, F (un ) → F (F (−1) (q)). Donc F (F (−1) (q)) ≥ q.
2. • Si F (x) ≥ q alors par définition F (−1) (q) ≤ x.
• Si x ≥ F (−1) (q) alors par croissance de F on a F (x) ≥ F (F (−1) (q)) donc
F (x) ≥ q par l’item 1.
31
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
p4
●
●
●
●
●
●
●
p3
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
p2
●
●
●
●
●
●
●
●
●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
p1
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
où on a utilisé l’item 2 pour la 2ème ligne, le fait que F est continue pour la
4ème ligne, et l’item 4 pour la dernière ligne. Comme ] − ∞, x] = ∩t>x ] − ∞, t[
on a P(F (Z) ≤ x) = limt→x,t>x P(F (Z) < t) = limt→x,t>x t = x. Donc F (Z) ∼
U [0, 1].
6. Soit q1 , q2 ∈ [0, 1] avec q1 ≤ q2 . Alors {x ∈ R : F (x) ≥ q2 } ⊂ {x ∈ R : F (x) ≥
q1 } donc F (−1) (q1 ) ≤ F (−1) (q2 ).
32
Remarque 2.8. Les item 1 et 4 peuvent "se déduire" à partir d’un dessin. On "voit"
également que les paliers de F correspondent à un point de discontinuité de F (−1) et
qu’un saut de F correspond à un palier de F (−1) .
Remarque 2.9. Dans un certain nombre de cas (cf exemples ci-dessous), la p-valeur
p(X) d’un test suit une loi uniforme sous H0 .
iid
Exemple 2.10. Un échantillon X1 , . . . , Xn ∼ N (µ, σ 2 ) avec σ inconnu. Problème
de test H0 : µ = µ0 contre H1 : µ ≤ µ0 . On utilise le test de Student φ =
1T ≤qT (n−1) . Alors d’après le théorème de Wassermann, la p-valeur observée s’écrit
α
p(xn1 ) = FT (n−1) (T (xn1 )). Donc la p-valeur p(X1n ), en tant que variable aléatoire, vé-
rifie p(X1n ) = FT (n−1) (T (X1n )). Or T (X1n ) ∼ T (n − 1) sous H0 . Donc, d’après l’item
5 de la proposition précédente, p(X1n ) suit une loi uniforme sous H0 .
n α o
p(xn1 ) = inf α ∈]0, 1[: |T (xn1 )| ≥ FT−1
(n−1) (1 − )
n
2
αo
= inf α ∈]0, 1[: FT (n−1) (|T (xn1 )|) ≥ 1 −
n h 2 io
n
= inf α ∈]0, 1[: α ≥ 2 1 − FT (n−1) (|T (x1 )|
= 2[1 − FT (n−1) (|T (xn1 )|].
Ainsi h i
p(X1n ) = 2 1 − FT (n−1) (|T (X1n )|)
Pour simplifier, on note T (X1n ) = T et FT (n−1) = F . On a alors si x ∈ [0, 1],
P 2[1 − F (|T |)] ≤ x
x
= P F (|T |) ≥ 1 −
2
x x
= P F (T ) ≥ 1 − , T ≥ 0 + P F (−T ) ≥ 1 − , −T ≥ 0
2
2
x
= 2P F (T ) ≥ 1 − , T ≥ 0
2
x
= 2P F (T ) ≥ 1 −
2
=x
On a utilisé
— pour la 3ème ligne : la symétrie de la loi de T .
33
— pour la 4ème ligne : 1− x2 > 1/2 si x ∈]0, 1[. Or F est la cdf de T (n−1), donc F
est continue et correspond à une loi symétrique. Donc F (x) = 1−F (−x). Donc
F (0) = 1/2 et comme F est strictement croissante, on a F (x) ≥ 12 =⇒ x ≥ 0
x
— pour la dernière ligne : 1 − F (T ) ∼ U [0, 1] et 1 − 2
∈]0, 1[.
Donc, à nouveau, p(X1n ) ∼ U [0, 1].
P (X ≤ qβ ) ≥ β et P (X ≥ qβ ) ≥ 1 − β.
Démonstration. 1. évident.
2. F (qβ ) ≥ β est l’item 1 de la proposition 2.7.
3. F (x− ) ≡ limt→x,t<x F (t) = limt→x,t<x P(] − ∞, t]) = P(] − ∞, x[). De plus si
x < qβ alors, par définition de qβ , on a F (x) < β. Donc F (qb− ) ≤ β.
P (X ≤ qβ ) ≥ β et P (X ≥ qβ ) ≥ 1 − β. (2.1)
34
Proposition 2.18. Soit X une variable aléatoire réelle de cdf F , et α ∈]0, 1[. Le
F
plus petit réel c tel que P(X > c) ≤ α est égal à q1−α .
On utilise une procédure de test φα = 1Pn Xi >c avec c choisi de façon à ce que le
i=1
niveau du test soit plus petit que α et tel que la région de rejet soit maximisée. On
B(n,1/2)
choisit donc c = q1−α .
Si on veut tester au niveau α
Théorème 2.20. Soit (Fn )n≥0 une suite de fonctions de répartition sur R et F une
fonction de répartition sur R. Alors Fn converge vers F en tout point de continuité
de F si et seulement si Fn(−1) converge vers F (−1) en tout point de continuité de
F (−1) .
Exemple 2.21. La loi de Student à n degrés de liberté tend vers la loi normale
standard. Φ−1 est continue. Donc, pour tout α ∈]0, 1[, qαT (n) → qαN (0,1) .
n→∞
On a besoin des quantiles pour les procédures de tests ainsi que pour les régions
de confiance. Parfois on ne sait pas calculer les quantiles de la loi mais on sait
simuler cette loi. Le quantile empirique peut alors être utilisé en remplacement du
vrai quantile.
On rappelle la notation suivante pour les statistiques d’ordre :
X(1) ≤ . . . ≤ X(n)
35
Intuition 2.2.1. Il s’agit donc des quantiles des cdf (lois) empiriques.
Proposition 2.23.
F̂n(−1) (β) = X(dnβe)
où on a noté dte = min{m ∈ N : m ≥ t}.
36
Théorème 2.27. Si F est dérivable en qβ avec F 0 (qβ ) > 0 alors
√ Loi β(1 − β)
n(q̂n,β − qβ ) −→ N 0, 2
F 0 (qβ )
H0 : F = F0 contre H1 : F 6= F0
hn (X1n , F0 ) = kF̂n − F0 k∞
Remarque 2.29. Il s’agit bien d’une statistique, c’est-à-dire que hn est bien mesu-
rable. En effet on peut montrer (grâce à la continuité à droite) que
hn (X1n , F0 ) = sup F̂n (x) − F0 (x).
x∈Q
hn (X1n , F0 ) = sup Mj
0≤j≤n
où, pour 1 ≤ j ≤ n − 1,
Mj = sup F̂n (x) − F0 (x)
x∈[X(j) ,X(j+1) [
37
et
M0 = sup F̂n (x) − F0 (x) et Mn = sup F̂n (x) − F0 (x)
x<X(1) x≥X(n)
et
−
M0 = sup 0 − F0 (x) = sup F0 (x) = F0 (X(1) )
x<X(1) x<X(1)
Et par la continuité de F0 ,
M0 = F0 (X(1) )
Considérons maintenant Mj pour 1 ≤ j ≤ n − 1. On a
j
Mj = sup − F0 (x).
x∈[X(j) ,X(j+1) [ n
j j
Mj = max{F0 (X(j+1) ) − , − F0 (X(j) )}.
n n
En rassemblant tous les résultats on obtient finalement
j j
hn (X1n , F0 ) = max max {F0 (X(j+1) )− }, max { −F0 (X(j) )}, F0 (X(1) ), 1−F0 (X(n) )
1≤j≤n−1 n 1≤j≤n−1 n
On obtient le résultat final en remarquant que
j j
max max { − F0 (X(j) )}, 1 − F0 (X(n) ) = max { − F0 (X(j) )}
1≤j≤n−1 n 1≤j≤n n
et
j j−1
max max {F0 (X(j+1) ) − }, F0 (X(1) ) = max max {F0 (X(j) ) − }, F0 (X(1) )
1≤j≤n−1 n 2≤j≤n n
j−1
= max {F0 (X(j) ) − }
1≤j≤n n
38
Définition 2.31. — On dit qu’une variable Z est diffuse si sa cdf est continue.
— "La statistique hn (X1n , F0 ) est libre de F0 " signifie que sa loi ne dépend pas de
F0 .
Proposition 2.34. Sous H0 , si F0 est continue alors hn (X1n , F0 ) est une statistique
libre de F0 et de loi continue.
Démonstration. Soit U1n = (U1 , . . . , Un ) est un n-échantillon iid de loi uniforme sur
iid (−1) (−1)
]0, 1[ Sous H0 , comme Xi ∼ F0 , on a (X1 , . . . , Xn ) ∼ (F0 (U1 ), . . . , F0 (Un ))
d’après la proposition 2.7. On a donc aussi, sous H0 ,
1 X
n
hn (X1n , F0 ) ∼ sup 1F (−1) (Ui )≤x − F0 (x)
x∈R n i=1 0
39
où on a noté G la fonction de répartition de la loi uniforme sur [0, 1], i.e. la fonction
définie par G(s) = s pour s ∈ [0, 1]. Cela montre que la loi, sous H0 , de hn (X1n , F0 )
est libre de F0 .
On prouve maintenant que la loi de hn (U1n , G) est continue. D’après la proposi-
tion 2.30, comme G est continue,
j j−1
hn (U1n , G) = max max{ − U(j) , U(j) − }
1≤j≤n n n
Comme la loi des Uj est absolument continue, celle de U(j) aussi (fait en TD1, on
a même ici U(j) ∼ Beta(j, n − j + 1), toujours d’après le TD1). Donc, d’après la
remarque 2.33, hn (U1n , G) est bien de loi continue.
iid iid
Exemple 2.35. Soit X1 , . . . , Xn ∼ N (0, 1) et Y1 . . . , Yn ∼ exp(1) alors
1 X
n 1 X
n 1X n
sup 1Xi ≤x − Φ(x) ∼ sup 1Yi ≤x − (1 − exp(−x)) ∼ sup 1Ui ≤s − s
x∈R n i=1 x>0 n i=1 s∈]0,1[ n i=1
Pour tout α ∈]0, 1[, si on note ξn,α le quantile d’ordre α de la loi de la statistique
hn (U1n , G), on a donc, par continuité de cette loi,
B(n, α) = { fonctions de répartitions G : ∀x ∈ R F̂n (x) − ξn,1−α ≤ G(x) ≤ F̂n (x) + ξn,1−α }
= {G : hn (X1n , G) ≤ ξn,1−α }
(En effet, la continuité de F0 n’est pas nécessaire pour obtenir cette égalité en loi).
Et comme
1 Xn 1 Xn
sup 1Ui ≤s − s ≤ sup 1Ui ≤s − s,
s∈Im(F0 ) n i=1 s∈]0,1[ n i=1
on a
P iid hn (X1n , F0 ) ≥ ξn,1−α ≤ P hn (U1n , G) ≥ ξn,1−α = α.
X1 ,...,Xn ∼ F0
40
Remarque 2.38. Quand le nombre de données n est grand, on utilise un test asymp-
totique.
On a aussi l’inégalité de Dvoretzky-Kiefer-Wolfowitz, et qui est valable sans
iid
condition sur F0 . Sous H0 : X1 , . . . Xn ∼ F0 ,
2
P(kF̂n − F0 k∞ > ) ≤ 2e−2n pour tout n ∈ N et tout > 0
41
On est alors ramené à prouver le résultat pour des variables uniformes sur [0, 1],
pour lesquelles les deux premiers problèmes ne se posent pas, puisque la cdf est ici
G(x) = x, définie sur le segment [0, 1], et continue. Il reste donc à prouver que,
presque sûrement, Ĝn converge uniformément vers G. Pour cela, il reste à régler le
dernier problème. Il suffirait, pour conclure à l’aide de Dini, de prouver qu’il existe
un ensemble mesurable A de probabilité 1 tel que, si ω ∈ A, alors pour tout x,
Ĝn (ω, x) tend vers G(x). On pose A = ∩q∈Q∩[0,1] A(q). Si ω ∈ A on a donc
→ G(q),
Ĝn (ω, q) n→∞ ∀q ∈ Q ∩ [0, 1]
De plus, comme Q est dénombrable, P(A) = 1. Il reste à prouver que, pour tout
s ∈ [0, 1], on a aussi Ĝn (ω, s) → G(s), si ω ∈ A. Soit donc ω ∈ A, s ∈ [0, 1] et
n→∞
> 0. Par densité de Q dans R, Il existe des rationnels q1 et q2 tels que s − ≤ q1 ≤
s ≤ q2 ≤ s + . Par croissance de Ĝn on a
Ĝn (q1 ) ≤ Ĝn (s) ≤ Ĝn (q2 )
Donc en passant à la limite sup et la limite inf (attention à ce stade on ne sait
pas encore que Ĝn (ω, s) converge, donc on doit utiliser les limites supérieures et
inférieures qui, elles, existent toujours), on obtient
inf Ĝn (ω, s) ≤ lim sup Ĝn (ω, s) ≤ G(q2 ) = q2 ≤ s + .
s − ≤ q1 = G(q1 ) ≤ lim n→∞
n→∞
Ces inégalités étant vraies pour tout > 0, on a bien limn→∞ Ĝn (ω, s) = s et la
preuve est terminée.
1 X
n
t
= sup 1Yi ≤t − (1 − e− Ȳ )
t>0 n i=1
42
P
t
La statistique supt>0 n1 ni=1 1Yi ≤t − (1 − e− Ȳ ) a une loi indépendante de λ.
On admet la continuité de cette loi.
Proposition 2.42. Sous H0 : F = G et si F est continue alors la loi de hn,m (X1n , Y1m )
ne dépend pas de F .
iid iid
Démonstration. Sous H0 , X1 , . . . , Xn , Y1 , . . . , Ym ∼ F donc, si U1 , . . . , Un , V1 , . . . , Vm ∼
U [0, 1], on a
F (−1) (U1 ), . . . , F (1 ) (Un ), F (−1) (V1 ), . . . , F (−1) (Vm ) ∼ (X1 , . . . , Xn , Y1 , . . . , Ym )
Ainsi on obtient, sous H0 ,
1 X
n
1 Xm
hn,m (X1n , Y1m ) = sup 1Xi ≤t − 1Yi ≤t
t∈R n i=1 m i=1
1 X
n
1 Xm
∼ sup 1F (−1) (Ui )≤t − 1F (−1) (Vi )≤t
t∈R n i=1 m i=1
1 X
n
1 Xm
= sup 1Ui ≤F (t) − 1Vi ≤F (t)
t∈R n i=1 m i=1
1 X
n
1 Xm
= sup 1Ui ≤s − 1Vi ≤s
s∈Im(F ) n i=1 m i=1
1 Xn
1 Xm
= sup 1Vi ≤s
p.s.
1Ui ≤s −
s∈]0,1[ n i=1 m i=1
43
On a utilisé la proposition 2.7 ainsi que la continuité de F . En effet, si F est conti-
nue, ]0, 1[⊂ Im(F ) ⊂ [0, 1] et on vérifie immédiatement que, presque sûrement, la
P P
fonction s 7→ n1 ni=1 1Ui ≤s − m1 mi=1 1Vi ≤s vaut 0 en s = 0 et s = 1.
où xn,m,1−α est le quantile d’ordre 1 − α de la loi de la statistique hn,m (X1n , Y1m ) sous
H0 .
Avec R : Pour tester si deux échantillons x et y ont la même loi, on peut utiliser
ks.test du package stat. La formule est ks.test(x,y).
44
aequo (tiré des documents pédagogiques de F-G Carpentier, cf biblio). L’échantillon
se nomme x.
> x= c(8.43, 8.70, 11.27, 12.92, 13.05, 13.05, 13.17, 13.44, 13.89,
18.90)
> ks.test(x,"pnorm",mean=13, sd=3)
Remarque 2.45. Pour les étudiants qui préfèrent Python à R, on peut appeler les commandes R
depuis Python. Par exemple on peut faire
Ensuite on utilise normalement la fonction qu’on a appelée rks, en prenant garde de transformer
aussi les entrées. Par exemple si on a un échantillon dans le vecteur x, et si on veut vérifier qu’il
s’agit d’un échantillon gaussien standard :
45
y=robjects.FloatVector(x)
z=rks(y,"pnorm")
On peut aussi utiliser directement les fonctions natives stats.kstest, qui se comportent différem-
ment en cas d’ex aequo.
46
2.5 Exercices
iid
Exercice 2.1. Soit X1 , . . . , Xn ∼ Be(p). On veut tester
est libre de F0 .
2. Proposer une procédure de test de
alors
β(F ) → 1
n→∞
H0 : F ∈ FN contre H1 : F 6∈ FN ,
où
FN = G : ∃ (µ, σ 2 ) ∈ (R × R∗+ ) tel que G = Nµ,σ2 .
47
Figure 2.2 – Exercice 3 : fonction de répartition empirique de l’échantillon (en
escalier) et fonction de répartition à tester.
48
1. Montrer que, sous H0 , Tn tend en loi vers la loi normale standard.
2. Montrer que l’erreur de première espèce du test de Student appliqué à l’échantillon (X1 , . . . , Xn )
tend vers α quand n tend vers l’infini.
Pour cela, on admettra le résultat suivant (qui est une généralisation du 2ème théorème de
Dini) : Si (Fn )n≥0 et F des fonctions de répartition, si F est continue et si Fn converge
simplement vers F alors la convergence est uniforme.
3. Montrer que la puissance du test tend simplement vers 1 quand n tend vers l’infini.
Exercice 2.6. L’objectif de cet exercice est de proposer une procédure de tests multiples lorsque
le nombre d’hypothèses à tester est élevé. On considère dans tout l’exercice (Ω, A, Pθ , θ ∈ Θ) un
modèle statistique.
Partie A. On se place tout d’abord dans le cadre simple où on veut tester
H0 : θ = θ0 contre H1 : θ ∈ Θ1
/ Θ1 . Pour cela, on dispose d’une observation réelle X de loi Pθ . Pour α ∈]0, 1[ donné , on
où θ0 ∈
considère un test de H0 contre H1 de la forme φα (X) = 1X≥kα où kα ∈ R. On note Fθ la cdf de
X sous Pθ . On suppose que Fθ est continue.
1. Montrer que la p-valeur observée de ce test s’écrit pour tout x ∈ R :
Θ0i ∩ Θ1i = ∅
On suppose pour simplifier que les hypothèses nulles sont des singletons Θ0i = {θ0i }. On note I0
l’ensemble des indices i pour lesquels H0i est vraie :
On cherche à construire une procédure de tests multiples qui retourne un ensemble R̂ ⊂ {1, . . . , m}
correspondant aux indices i pour lesquels H0i est rejetée. On note FP le cardinal de l’ensemble
des indices correspondant aux hypothèses nulles rejetées à tort et TP le cardinal de l’ensemble des
indices correspondant aux hypothèses nulles rejetées à raison :
FP = card(R̂ ∩ I0 ), TP = card(R̂ \ I0 ).
FP est le cardinal des faux positifs et TP celui des vrais positifs. Idéalement, on cherche une
procédure de tests de sorte que FP soit petit et TP soit grand. On note p̂i la p-valeur du test de
H0i contre H1i . Donc p̂i est une statistique satisfaisant, pour tout u ∈]0, 1[,
49
(b) En utilisant (2.6), en déduire que
P(FP > 0) ≤ α.
2. La procédure de Bonferroni contrôle le nombre de faux positifs mais peut produire un trop
petit nombre de vrais positifs. On dit que c’est une procédure trop conservative. Aussi, on
propose l’alternative suivante. On se donne une fonction f : {0, 1, . . . , m} → [0, m] supposée
croissante et on ordonne les statistiques p̂i par ordre croissant :
avec
αf (k)
k̂ = max k ∈ {1, . . . , m} : p̂(k) ≤ .
m
αf (k)
En particulier, on pose k̂ = 0 et R̂ = ∅ si pour tout entier k, p̂(k) > m .
(a) Montrer que
k̂ = card(R̂)
et que pour j ≥ k̂,
f (k̂) ≤ f (min(j, m)).
(b) Établir alors que
X 1{k̂≥1}
FDR = E 1{p̂i ≤αf (k̂)/m} × .
i∈I0
k̂
alors
FDR ≤ α.
(f) Donner un exemple de fonction f satisfaisant (2.7).
50
Remarque : on peut aussi généraliser ces résultats au cas d’hypothèses nulles composites (cf
examen 2014 ou partiel 2018).
Remarque : si les m tests sont indépendants, on peut en fait prendre f égale à l’identité : alors
on a plus de vrais positifs que dans le cas précédent tout en ayant quand même un FDR borné par
α. La procédure est alors la suivante :
( )
αk̂
R̂ = i ∈ {1, . . . , m} : p̂i ≤
m
avec
αk
k̂ = max k ∈ {1, . . . , m} : p̂(k) ≤ .
m
Elle s’appelle la procédure de Benjamini-Hochberg (cf partiel 2018).
Code R : si on a calculé les p-valeurs des m tests indépendants dans le vecteur p alors on
peut utiliser le code suivant pour calculer R̂, l’ensemble des indices des hypothèses rejetées par la
procédure de Benjamini-Hochberg quand on, veut un FDR plus petit que 5% :
k<-sum(sort(p)<=0.05*(1:m)/m) # k chapeau
R<-(1:m)[p<=0.05*k/m]#
Il existe un certain nombre de méthodes basées sur les p-valeurs pour résoudre ce type de
problème de tests multiples. Par exemple on peut citer la procédure de Berk-Jones modifiée.
51
Chapitre 3
Tests robustes
3.1 Un exemple
Un exemple de question à laquelle on souhaite répondre dans ce chapitre est
la suivante : les hommes gagnent-ils plus que les femmes ? Pour répondre à cette
question, imaginons que nous disposions d’un échantillon X1n1 = (X1 , . . . , Xn1 ) de
salaires de femmes et d’un échantillon Y1n2 = (Y1 , . . . , Yn2 ) de salaires d’hommes.
Nous ferons des tests différents selon que
1. Les échantillons sont iid et indépendants entre eux i.e.
iid iid
X1 , . . . , X n 1 ∼ X et Y1 , . . . , Yn2 ∼ Y X1n1 ⊥⊥ Y1n2
52
test sur un paramètre de position. Deux exemples usuels de paramètres de position
sont la moyenne et la médiane. On pourrait traduire le fait que les femmes gagnent
autant que les hommes par "la variable différence Y1 − X1 a une moyenne égale à 0",
ou bien, si on préfère utiliser la médiane, on pourrait le traduire par "la médiane de
la différence est égale à 0".
Autrement dit, si nous choisissons la moyenne comme paramètre de position,
H0 : la moyenne de Y1 − X1 est égale à 0
contre
H1 : la moyenne de Y1 − X1 est strictement positive
Problèmes éventuels qu’on peut avoir pour réaliser ce test dans la pratique :
— l’échantillon n’est pas de loi normale,
— les variables sont gaussiennes mais pas de même variance : par exemple on
peut avoir Xi ∼ N (µ, σi2 ).
— l’échantillon est contaminé par des outliers (=observations aberrantes)
53
Code R
Le test de Student sur un échantillon dans R peut se faire par la procédure
t.test.
Le test de Shapiro-Wilk peut se faire avec shapiro.test.
et on veut tester :
H0 : µ1 = µ2 contre H1 : µ1 6= µ2
(ou bien H1 : µ1 < µ2 ou bien H1 : µ1 > µ2 )
On note σ 2 la variance commune et on suppose que σ est inconnu.
On utilise alors la variable
V̄ − Ū
T = q1 1
σ̂ n + p
où on a posé
Xn Xp
1
σ̂ 2 = (Ui − Ū )2 + (Vi − V̄ )2 .
n + p − 2 i=1 i=1
V̄ − Ū
q
1 1
∼ N (0, 1)
σ n
+ p
et
σ̂ 2 χ2 (n + p − 2)
∼ .
σ2 n+p−2
De plus
n
X p
X
2
(Ui − Ū ) ⊥⊥ Ū , (Vi − V̄ )2 ⊥⊥ V̄ , U1n ⊥⊥ V1p .
i=1 i=1
Donc
σ̂ 2 ⊥⊥ Ū − V̄
54
Et finalement
V̄ −Ū
q
1
σ n
+ p1
T = σ̂ ∼ T (n + p − 2).
σ
55
par exemple deux individus de même âge. Il peut aussi s’agir du même individu, à
qui on a donné deux traitements différents à deux moments différents.
De manière générale, quand on considère des échantillons appariés, cela signifie
que
— soit Ui et Vi correspondent à une mesure sur le même individu,
— soit les individus sont différents mais ils sont regroupés en fonction de cova-
riables (sexe, âge etc).
Les tests pour données appariées sont essentiellement basés sur le fait de prendre
la différence des deux mesures et ensuite de faire un test sur l’ échantillon résultant.
Xi = Ui − Vi , i = 1, . . . , n
On suppose que
H0 : µ = 0, contre H1 : µ 6= 0
où
√ Ū − V̄ n
1 X
T = n σ̂ 2 = (Ui − Vi − Ū − V̄ )2
σ̂ n − 1 i=1
.
Une des hypothèses faites est que la distribution des Xi est la même pour tout i. En particulier
la variance doit être la même pour tout i. Certains auteurs recommandent de faire une vérification
graphique de cela avec un graphe de "Bland-Altman". Fréquemment la dispersion est proportion-
nelle au niveau et une transformation logarithmique est utile pour remédier à ce problème.
Avec R
On peut utiliser t.test et il faut préciser paired=T pour dire que les données
sont appariées : t.test(x,y,paired=T,var.equal=T).
De façon équivalente on peut utiliser t.test(x-y,var.equal=T).
56
3.2.4 Importance des conditions d’application
Le test de Student
A retenir : les tests non paramétriques s’appliquent dans des situations plus gé-
nérales et sont donc plus robustes. On les utilise en général quand les conditions
d’application des tests paramétriques ne sont pas vérifiées (ou pas vérifiables). Tou-
tefois, un test paramétrique peut devenir performant avec une grande taille d’échan-
tillon même si les conditions théoriques d’application du test paramétrique ne sont
pas exactement vérifiées. En particulier, pour le test de Student de comparaison
de deux populations, quand les tailles des échantillons sont importantes et sous des
conditions assez faibles sur la loi des échantillons, le test de Student est valide, même
si les échantillons ne sont pas gaussiens. Ce résultat est à rapprocher de ce qui se
produit en modèle linéaire quand les erreurs ne sont pas gaussiennes.
De manière générale, on peut cependant préférer utiliser systématiquement les
tests de Wilcoxon quand on ne sait pas si échantillon sont gaussiens. En effet, la
performance des tests de Wilcoxon effectué sur des échantillons gaussiens n’est pas
tellement moins bonne que la performance des tests de Student. De plus, les tests
de Wilcoxon ont souvent une meilleure performance que celle du test de Student
quand l’échantillon n’est pas gaussien (même avec une grande taille d’échantillon). .
Le reste de la section 3.2.4 est facultatif
Revenons à notre exemple lié aux salaires. Supposons ici que les échantillons sont iid et in-
dépendants entre eux. Nous avons donc un échantillon iid de salaires de femmes, de taille n1 ,
et un échantillon iid de salaires d’hommes, de taille n2 , et ces deux échantillons sont supposés
indépendants. Supposons que les échelles des distributions de salaires soient les mêmes (même
dispersion).
Si les données normales alors nous choisirons le test de Student. Que se passe-t-il si nous nous
trompons et que nous appliquons le test de Student à deux échantillons de loi non gaussienne par
exemple ? Ou s’il y a des outliers dans les échantillons ?
Dans les simulations qui suivent, nous utilisons aussi le test de Wilcoxon de la somme des
rangs pour comparer deux échantillons indépendants (appelé aussi "test de Mann-Whitney"). Ce
test non paramétrique peut aussi être utilisé pour comparer les positions de deux populations.
Comme tout test non paramétrique, il a des conditions d’application beaucoup plus générales
que les tests paramétriques. En particulier pour l’appliquer, il n’est pas nécessaire d’avoir des
échantillons gaussiens. Ce test sera étudié dans la suite.
Nous voulons donc illustrer ici ce qui se produit quand on n’est pas dans les conditions d’ap-
plication du test de Student (et ensuite comparer sa performance avec le test de Mann-Whitney)
pour montrer l’intérêt des tests non paramétriques. Plus précisément nous allons simuler des va-
riables de lois normales et aussi des lois non normales : nous regardons ce qui se passe pour des
échantillons de loi de student à 3 degrés de liberté et de loi de Cauchy . Comme la loi de Cauchy
n’a pas de moyenne, ce que nous utilisons comme paramètre de position pour la loi de Cauchy est
sa médiane.
57
sim=function(type,n,a)
# fonction qui simule un échantillon de taille n
{
return(switch(type,norm=rnorm(n,mean=a),cauchy=rcauchy(n,a), student3=rt(n,df=3,a)))
}
test=function(n,type,a,b,n1,n2,outliers=F)
#simulation, calcul de la p-valeur de chaque test
{
u=rep(0,n); v=rep(0,n); w=rep(0,n);
#fait avec lapply par principe mais boucle for pas plus lente ici
lapply(1:n,function(i)
# on fait n simulations et les 3 tests sur chaque simulation
{
x=sim(type,n1,a); y=sim(type,n2,b);# simulation des deux échantillons
if (outliers) {x[1:10]=rnorm(10,3)}
#calcul de la p-valeur de chaque test
v[i]<<-t.test(x,y,var.equal=T)$p.value; # Student
w[i]<<-wilcox.test(x,y)$p.value; # Mann-Whitney
})
# on regarde le taux d'erreurs si on est sous H0 et la puissance si on est sous H1
return(list( "Student"=sum(v<0.05)/n,
"Mann-Whitney"=sum(w<0.05)/n))
## $Student
## [1] 0.0513
##
## $`Mann-Whitney`
## [1] 0.0503
#sous H1
test(10000,"norm",1,2,30,30)
## $Student
## [1] 0.9686
##
## $`Mann-Whitney`
## [1] 0.9613
#loi de Student à 3 degrés de liberté (même variance)
#sous HO
test(10000,"student3",1,1,50,50)
## $Student
## [1] 0.0443
46
##
## $`Mann-Whitney`
## [1] 0.0466
#sous H1
test(10000,"student3",1,2,50,50)
## $Student
## [1] 0.8839
##
## $`Mann-Whitney`
## [1] 0.979
#loi de Cauchy, petits échantillons
#sous H0
test(10000,"cauchy",1,1,15,15)
## $Student
## [1] 0.0218
##
##$`Mann-Whitney`
##[1] 0.0446
#sous H1
test(10000,"cauchy",1,2,15,15)
## $Student
## [1] 0.0728
##
## $`Mann-Whitney`
## [1] 0.2701
#loi de Cauchy, grands échantillons
#sous H0
test(10000,"cauchy",1,1,100,100)
## $Student
## [1] 0.0201
##
## $`Mann-Whitney`
## [1] 0.0495
#sous H1
test(10000,"cauchy",1,2,100,100)
## $Student
## [1] 0.0762
##
## $`Mann-Whitney`
## [1] 0.9548
# lois normales avec 10\% d'outliers "assez gros"
# sous H0
test(10000,"norm",1,1,100,100,outliers=T)
## $Student
## [1] 0.2278
##
## $`Mann-Whitney`
47
## [1] 0.1572
#sous H1
test(10000,"norm",1,2,100,100,outliers=T)
## $Student
## [1] 0.9998
##
## $`Mann-Whitney`
## [1] 0.9998
#la même chose mais sans outliers
# sous H0
test(10000,"norm",1,1,100,100)
## $Student
## [1] 0.0456
##
## $`Mann-Whitney`
## [1] 0.0441
#sous H1
test(10000,"norm",1,2,100,100)
## $Student
## [1] 1
##
## $`Mann-Whitney`
## [1] 1
• Le test de Mann-Whitney est souvent moins performant que le test de Student quand on est dans les
conditions d’application du test de Student, mais la différence est souvent faible.
• Avec des lois de Student à 3 degrés de liberté, le test de Mann-Whitney est plus performant que le test
de Student. Le test de Student est valide dans le cas où les échantillons sont de taille suffisante. Ce
comportement peut grossièrement s’expliquer ainsi : la variance est finie, donc quand n1 et n2 sont
suffisamment grands, le test fonctionne assez bien (cf aussi TD2 exo 5) .
• Le test de Student se comporte mal avec des échantillons de loi de Cauchy, même si les tailles des
échantillons sont grandes. Ce comportement peut grossièrement s’expliquer ainsi : les queues de la loi
de Cauchy sont si lourdes que la variance est infinie (même la moyenne est infinie dans ce cas).
• La performance du test de Student est plus affectée par la présence d’outliers que le test de Mann-
Whitney.
Dans ces deux dernières situations, c’est-à-dire quand il y a présence d’un grand nombre d’outliers ou quand
loi est à queues lourdes, le test de Mann-Whitney que l’on va introduire dans la suite se comporte mieux que
le test de Student. On dit qu’il est plus robuste.
48
3.3 Test du signe
Nous venons de voir un test paramétrique, le test de Student, qui peut être utilisé
pour comparer deux populations indépendantes ou bien comparer deux traitements
sur des données appariées. Ce test repose sur le caractère gaussien des données.
On va construire maintenant des tests reposant sur des hypothèses beaucoup
plus faibles sur les données.
On commence par le test du signe et le test de Wilcoxon des rangs signés, qu’on
peut plus ou moins voir comme des versions non-paramétriques du test de Student.
∀x ∈ R, P(U = x) = 0
Cela revient à dire que sa distribution est continue, c’est-à-dire que sa cdf est conti-
nue
Les conditions :
1. les Xi sont indépendantes
2. Les Xi ont une médiane commune m.
3. P (Xi = m) = 0
Remarquez que les Xi ne sont pas nécessairement identiquement distribuées.
L’hypothèse nulle est :
H0 : m = 0
Remarquons que m = 0 implique ici que P(Xi ≤ 0) = 1/2.
En effet si 0 est la médiane commune des Xi alors (cf chapitre 2)
61
(c’est donc ici qu’intervient la condition 3. )
On pose
Yi = 1Xi ≤0 .
Faisons d’abord, pour simplifier l’exposition, l’hypothèse que les Xi sont de même
loi.
On a
iid
Yi ∼ Be(p), avec p = P(Xi ≤ 0).
Donc H0 se réécrit
H0 : p = 1/2.
Donc on se ramène à un test d’égalité sur le paramètre p d’un échantillon iid de
v.a. de Bernoulli Yi .
Imaginons que l’alternative soit la suivante :
H1 : m 6= 0
H1 : p 6= 1/2.
φα (Y1 , . . . , Yn ) = 1 P
n B(n,1/2)
| Y −n/2|>q1− α
i=1 i
−n
2
2
φα (X1 , . . . , Xn ) = 1 P
n B(n,1/2) n
| 1
i=1 Xi ≤0
−n/2|>q1− α −2
2
H1 : m < 0
φα (X1 , . . . , Xn ) = 1 P
n B(n,1/2
i=1
1{Xi ≤0} >q1−α
Pn Pn
Si l’hypothèse alternative est H1 : m > 0, on remplace i=1 1Xi ≤0 par i=1 1Xi ≥0
dans la formule ci-dessus.
Maintenant, que se passe-t-il si les Xi ne sont pas de même loi ? Pour simplifier,
supposons que H1 corresponde au fait que la médiane commune est strictement
négative :
H1 : m < 0
Alors "ça marche quand même " : en effet, sous H0 on a bien, du fait que 0 est la
médiane commune,
iid
Yi ∼ Be(1/2)
62
donc si on pose
φα (X1 , . . . , Xn ) = 1 P ,
n B(n,1/2)
i=1
1Xi ≤0 >q1−α
Avec R
On peut utiliser la procédure binom.test qui fait un test sur le paramètre p
d’un échantillon de va de Bernoulli. Si l’échantillon se trouve dans un vecteur x, et
si a une alternative bilatéral H1 : m 6= 0, on peut utiliser la commande suivante
binom.test(sum(x>0),n=length(x),p=0.5,alternative="two.sided")
Pour l’alternative H1 : m < 0 on met alternative="less"
Pour H1 : m > 0 on met alternative="greater".
63
Le test est donc
φα (U1 , . . . , Un , V1 , . . . , Vn ) = 1 P
n B(n,1/2)
i=1
1{Ui ≤Vi } >q1−α
Si au contraire, on pense que soit les deux médicaments ont la même efficacité,
soit le premier est plus efficace, alors on échange juste les rôle de Ui et Vi , ce qui
donne le test
φα (U1 , . . . , Un , V1 , . . . , Vn ) = 1 P
n B(n,1/2)
i=1
1{Vi ≤Ui } >q1−α
Si on n’a pas d’a priori sur les médicaments, c’est-à-dire si on ne sait pas quel
médicament est susceptible d’être plus efficace, l’alternative est alors H1 : m 6= 0 et
on fait le test
φα (X1 , . . . , Xn ) = 1 P
n B(n,1/2)
| 1
i=1 Ui ≤Vi
−n/2|>q1− α −n
2
2
Remarque 3.2. Le test du signe n’utilise que très peu d’information sur les variables
Ui − Vi (uniquement leur signe, pas leurs valeurs absolues). C’est donc un test peu
puissant. Quel est alors l’intérêt de parler du test du signe ? Il se peut que les signes
des Ui − Vi soit la seule donnée disponible : c’est en effet le cas si la question posée
aux patients qui ont testé les deux médicaments est "quel est le meilleur des deux ?"
(au lieu de noter les médicaments sur une échelle de 1 à 10 par exemple).
Remarque 3.3. Concrètement, il faut bien vérifier que la valeur 0 n’est pas dans
l’échantillon.
Avec R
Comme il s’agit du test du signe sur les variables Xi , on utilise exactement la
même procédure que dans le cas d’un seul échantillon. Il suffit donc de calculer
l’échantillon des différences puis de faire le test sur cet échantillon à l’aide de la
fonction binom.test.
(Remarque 26 enlevée)
Remarque 3.4. On pourrait utiliser ce test pour tester que les échantillons sont
de même loi. Supposons que l’on nous donne deux échantillons indépendants U1n =
(U1 , . . . , Un ) et V1n = (V1 , . . . , Vn ). Supposons que les Ui sont iid de loi F continue,
et les Vi sont iid de loi G continue. Nous voulons donc tester l’égalité des lois :
H0 : F = G contre H1 : F 6= G.
64
Pour cela nous voulons utiliser le test du signe. Alors nous devons vérifier si, sous
l’hypothèse F = G, l’hypothèse nulle associée au test du signe est vérifiée, c’est-à-
dire si m = 0 en notant m la médiane de V − U . Comme V − U est diffuse, cela
revient à montrer que
P(U ≤ V ) = P(V ≤ U )
Or on a, sous l’hypothèse F = G,
(U, V ) ∼ (V, U )
U −V ∼V −U
H0 : U ∼ V
alors une p-valeur observée grande n’implique pas que H0 est vraie. Par exemple, il
n’est pas rare d’avoir une p-valeur grande si on fait le test du signe sur un échantillon
de loi normale et l’autre de loi de Cauchy, alors que l’on est bien sous H1 : F 6= G.
Le test du signe, vu comme un test d’homogénéité, est donc un exemple de test
particulièrement peu puissant : "il ne voit" pas certaines alternatives.
65
3.4 Statistiques d’ordre et de rang
Définition 3.5. Soient X1 , . . . , Xn n v.a. réelles. La statistique d’ordre (X(1) , . . . , X(n) )
est définie par
{X(1) , . . . , X(n) } = {X1 , . . . , Xn }
et
X(1) ≤ X(2) ≤ . . . ≤ X(n)
On pose
X ∗ = (X(1) , . . . , X(n) )
66
Démonstration. Pour tout i 6= j on a
H0 : m = 0
67
Remarque 3.8. On a 0 ≤ Wn+ ≤ n(n+1) 2
. Le cas Wn+ = 0 correspond à tous les
Xi < 0, le cas Wn+ = n(n+1)
2
correspond au cas où tous les Xi > 0.
P
Si on pose en plus Wn = ni=1 R|X| (i)1{Xi <0} et si P(Xi = 0) = 0 (ce qui est le
−
n(n + 1)
Wn+ + Wn− =
2
Expliquons rapidement l’idée derrière ce test. Supposons pour fixer les idées que
H1 : m > 0.
L’idée est que, sous H1 , il y a plus de Xi positifs que de Xi négatifs. Jusque là c’est
même idée que pour le test du signe. Mais en plus, du fait de la symétrie, les Xi
positifs ont tendance à être plus grands en valeur absolue que les Xi négatifs, c’est
là qu’on utilise une information supplémentaire par rapport au test du signe. Donc
sous cette alternative, Wn+ sera "grand".
Evidemment si l’alternative est H1 : m < 0, alors Wn+ sera au contraire "petit".
68
●
3
● ●
●
2
2
● ●
●
●
● ●
● ●
●
●
1
1
● ● ●
●
●
● ● ● ●
●
●
x
y
● ● ●
● ● ●
●●
● ●
●
● ● ● ●●
0
0
● ●
● ●
●
● ● ● ●
●
−1
● ● −1
−2
−2
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Index Index
●
6
●
4
●
●
●
● ●
x
● ●
● ● ● ●
●
2
● ●
●
●
● ●
● ●
● ● ●
● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ●
● ● ●● ● ● ●● ● ●
●
●●●●
0
● ● ●●
● ● ● ● ● ● ●
● ● ● ● ● ●
● ● ● ●
● ● ●
● ● ● ● ●
●● ● ● ●
●● ● ●●
● ● ●
●
0 20 40 60 80 100
Index
69
5. Asymptotiquement, on a
On admet la preuve. Cependant pour les étudiants intéressés, voici une partie
de la preuve.
Démonstration. On se place sous H0 .
1. On a
n
X
Wn+ = R|X| (i)1Xi >0
i=1
Xn
= j1{Xσ >0}
|X| (j)
j=1
−1
où on a noté σ|X| = R|X| . De même
n
X
Wn− = j1{Xσ <0}
|X| (j)
j=1
X1n ∼ −X1n
f (X1n ) ∼ f (−X1n ).
70
et donc n
X
E(Wn+ ) = jP(Xσ|X| (j) > 0) = n(n + 1)/4
i=1
où
iid
Y1 , . . . , Yn ∼ Be(1/2).
En effet, soit (1 , . . . , n ) ∈ {0, 1}n , on a
P (1Xσ|X| (1) >0 , . . . , 1Xσ|X| (n) >0 ) = (1 , . . . , n )
X
= P (1Xs(1) >0 , . . . , 1Xs(n) >0 ) = (1 , . . . , n ), σ|X| = s
s∈Sn
X
= P (1Xs(1) >0 , . . . , 1Xs(n) >0 ) = (1 , . . . , n ) P(σ|X| = s)
s∈Sn
1 X
= P(σ|X| = s)
2n s∈Sn
1
= n
2
= P (Y1 , . . . , Yn ) = (1 , . . . , n )
On a utilisé
• ligne 2 : les probabilités totales.
• ligne 3 : l’indépendance entre (|X1 |, . . . , |Xn |) et (1X1 >0 , . . . , 1Xn >0 ) en-
traine l’indépendance entre (|X1 |, . . . , |Xn |) et (1Xs(1) >0 , . . . , 1Xs(n) >0 ) car
s est fixe ! (et σ|X| est une fonction de |X).
• ligne 4 : s est fixe et les variables X1 , . . . , Xn sont indépendantes donc les
variables Xs(1) , . . . , Xs(n) sont indépendantes. Donc
P 1Xs(1) >0 , . . . , 1Xs(n) >0 ) = (1 , . . . , n )
= P(1Xs(1) >0 = 1 ) . . . P(1Xs(n) >0 = n )
71
Pn
Ceci prouve l’item 3 : en effet, Wn+ ∼ j=1 jYj . Ceci permet aussi de trouver
la valeur de la variance : en effet
n
X n
X
Var( j1{Xσ (j) >0}
) = j 2 Var(Yj )
|X|
j=1 j=1
n
X j2 n(n + 1)(2n + 1)
= =
j=1 4 24
n(n + 1)
Wn+ ∼ − Wn+
2
c’est-à-dire
Wn+ ∼ 2b − Wn+
avec
n(n + 1)
b=
4
En conséquence, le test, pour l’alternative H1 : m 6= 0, est
α
où q1− α2 est le quantile d’ordre 1 − 2
de la loi de Wn+ sous H0 (il s’agit de la loi de
Pn iid
j=1 jYj avec Y1 , . . . , Yn ∼ Be(1/2) d’après la preuve).
Avec R
Pour tester H0 : m = 0 contre H1 : m 6= 0, si l’échantillon se trouve dans
un vecteur x, on peut utiliser wilcox.test(x,alternative="two.sided"). Pour
H1 : m > 0, on met alternative="greater" et pour H1 : m < 0, on met
alternative="less".
72
3.5.2 Echantillons appariées
On suppose disposer de deux échantillons appariés de taille n : (U1 , . . . , Un ) et
(V1 , . . . , Vn ) . On veut savoir si l’un des échantillons a "tendance à prendre des valeurs
plus grandes que l’autre" (penser à l’exemple des traitements médicamenteux). C’est
la même problématique que pour le test du signe sur deux échantillons. On modélise
le problème de la même manière.
Comme pour le test du signe, on pose
Xi = Vi − Ui
et nous utilisons le test de Wilcoxon des rangs signés sur l’échantillon X1n .
Nous supposons donc que
— Les Xi sont indépendants entre eux (mais pas forcément de même loi).
— Les Xi sont diffuses.
— Les Xi ont une médiane commune m.
— Les Xi sont de loi symétrique par rapport à m.
L’hypothèse nulle est
H0 : m = 0
Si on pense que soit on est sous H0 , soit les Ui prennent des valeurs plus petites
que les Vi , alors on pose comme hypothèse alternative
H1 : m > 0
Remarque 3.14. Comme on l’a déjà remarqué pour le cas du test sur la médiane
d’un seul échantillon, le test des rangs signés est plus puissant que le test du signe.
Donc il est conseillé d’utiliser le test des rangs signés plutôt que le test du signe si
on peut le faire, c’est-à-dire si on a accès aux valeurs de l’échantillon des Xi et pas
seulement à leur signe et si la loi des Xi est symétrique.
Avec R
On peut utiliser la fonction wilcox.test. Si nos échantillons sont dans des vec-
teurs x et y, et si H1 : m 6= 0, on écrit
wilcox.test(x,y,paired=T, alternative="two.sided").
ou son équivalent wilcox.test(x-y, alternative="two.sided").
(Mêmes changements possibles d’alternative que précédemment. )
On peut aussi avoir des données correspondants à deux colonnes d’un dataframe.
Un exemple : on veut savoir si les salaires des hommes d’une entreprise sont du même
ordre que les salaires des femmes ou bien plus élevés. On suppose que l’on a apparié
les données (par exemple on a rassemblés les salaires selon l’âge de la personne).
On suppose que ces salaires apparaissent dans un dataframe nommé salaires avec
pour colonnes femmes et hommes, on peut alors utiliser
73
Remarque 3.15. Il faut alors faire attention au problème des ex aequo ("ties" en
anglais) quand on utilise la procédure wilcox.test. On peut quand même faire le
test, mais il n’est jamais exact, c’est-à-dire qu’il repose automatiquement sur une
approximation gaussienne.
Comme l’approximation gaussienne n’est valable que pour des grands échan-
tillons, on ne peut pas trop se fier au résultat de wilcox.test quand il y a des
ex aequo et quand la taille d’échantillon est trop petite (pas de problème en revanche
avec les éventuels ex aequo si la taille est suffisamment grande).
Dans le cas d’ex aequos, on reçoit le message suivant warning message : cannot
compute exact p-values with ties.
Ce qu’on entend par ex aequo ici, c’est un ex aequo dans l’échantillon xn1 des
différences ou un ex aequo dans l’échantillon des valeurs absolues des différences
(|x1 |, . . . , |xn |).
Remarque 3.16. Pour tester l’hypothèse de symétrie des Xi quand il s’agit d’un échantillon i.i.d,
on peut par exemple commencer par représenter les données (histogramme par exemple ou densité
cf chap4), (ou utiliser un des nombreux tests de symétrie : par exemple le test symmetry.test du
package lawstat ( attention cette fonction est un peu lente). En cas d’asymétrie sur l’échantillon
des différences il semble préférable d’utiliser le test du signe.
Certains praticiens utilisent une transformation des Xi pour rendre l’échantillon symétrique (mais le test fait
sur l’échantillon transformé n’est alors pas un test sur la médiane des Xi )
P(X1 < X2 < X3 ) = P(X1 < X3 < X2 ) = P(X2 < X1 < X3 ) = P(X2 < X3 < X1 )
= P(X3 < X2 < X1 ) = P(X3 < X1 < X2 ) = 1/6
74
On montre maintenant que σ et X ∗ sont indépendantes. On veut montrer que,
pour tout borélien B de Rd et toute permutation s de Sn , on a
P(X ∗ ∈ B ∩ σ = s) = P(X ∗ ∈ B)P(σ = s).
Et comme P(σ = s) = n!1 , cela revient à montrer que, pour toute permutation
s ∈ Sn
P(X ∗ ∈ B) = n!P(X ∗ ∈ B ∩ σ = s)
Comme les Xi sont indépendantes et de même loi, elles sont interchangeables et
donc, pour tout s et tout B,
∗
P X ∈ B ∩ σ = s = P Xs(1) < . . . < X(s(n) , (Xs(1) , . . . , Xs(n) ) ∈ B
= P X1 < . . . < Xn , (X1 , . . . , Xn ) ∈ B)
D’autre part, le théorème des probabilités totales permet d’écrire :
X
P(X ∗ ∈ B) = P(X ∗ ∈ B ∩ σ = s)
s∈Sn
X
= P X1 < . . . < Xn , (X1 , . . . , Xn ) ∈ B)
s∈Sn
= n!P X1 < . . . < Xn , (X1 , . . . , Xn ) ∈ B)
= n!P X ∗ ∈ B ∩ σ = s
où s est une permutation quelconque.
La principale conséquence de ce théorème est que la loi de RX ne dépend pas de
la loi des Xi . On en déduit que toute variable aléatoire qui ne s’exprime qu’à l’aide
du vecteur de rangs d’observations i.i.d. de loi continue a une loi indépendante
de ces observations. C’est bien ce que l’on cherche à obtenir en statistique non
paramétrique, où la loi des observations n’appartient pas à une famille paramétrée
connue. On pourra donc faire de l’estimation et des tests non paramétriques à l’aide
des rangs des observations.
Remarque 3.18. Pour tout s fixé (non aléatoire) dans Sn on a
(Xs(1) , . . . , Xs(n) ) ∼ (X1 , . . . , Xn )
Mais ça n’est pas vrai si la permutation est aléatoire (à moins qu’elle ne soit indé-
pendante des Xi ). Par exemple, on a évidemment
X ∗ = (Xσ(1) , . . . , Xσ(n) ) (X1 , . . . , Xn )
Proposition 3.19. Soient X1 , . . . , Xn n v.a. i.i.d. de loi continue de vecteur des
rangs RX = (R1 , . . . , Rn ) . Et, pour tout entier s tel que 1 ≤ s ≤ n, et pour toute
suite d’entiers distincts (r1 , . . . , rs ) dans {1, . . . , n}, on a
1
P (R1 , . . . , Rs ) = (r1 , . . . , rs ) = .
n(n − 1) . . . (n − s + 1)
En particulier, pour tout i ∈ {1, . . . , n}, Ri suit une loi uniforme sur {1, . . . , n}.
75
Démonstration. Pour simplifier, considérons d’abord trois cas simples.
— s = n : P (R1 , . . . , Rn ) = (r1 , . . . , rn ) = P(R = r) où r = (r1 , . . . , rn ) ∈ Sn .
Donc, d’après le théorème précédent, on a
1
P (R1 , . . . , Rn ) = (r1 , . . . , rn ) = ,
n!
ce qui est bien le résultat annoncé.
— s = 1 : P(R1 = s) = P( « X1 est le s-ème plus petit élément de l’échantillon ») =
1
n
, toujours du fait que les Xi sont interchangeables.
— s = 2 : alors
1 1
P (R1 , R2 ) = (s1 , s2 ) = P(R1 = s1 )P(R2 = s2 | R1 = s1 ) =
nn−1
H0 : F = G
Ce n’est donc pas un test sur la médiane ou la moyenne contrairement aux tests
précédents ( signe et Wilcoxon des rangs signés). Cependant l’alternative n’est pas
F 6= G. D’ailleurs si vous faites le test avec la commande R associée, il sera écrit alternative
hypothesis: true location shift is not equal to 0.
On suppose en fait que U et V ont la même loi, à un paramètre de position près,
c’est-à-dire
— soit U et V ont la même loi (H0 )
— soit U a tendance à prendre des valeurs plus grandes que V , ou le contraire
(H1 ).
Autrement dit
76
Exemple 3.20. veut tester un nouveau médicament par rapport à un ancien mé-
dicament. On donne le premier à un groupe de n personnes, et le deuxième à un
groupe de p personnes, ces deux groupes étant cette fois-ci indépendants. On veut
voir si le nouveau médicament est plus efficace que l’ancien.
Remarque 3.21. On peut voir le test de Student comme un test d’égalité en loi,
quand on suppose que les données sont gaussiennes et ne peuvent différer (éventuel-
lement) que par leur moyenne. En ce sens le test de Mann-Whitney peut être vu
comme vu comme une version non-paramétrique (et plus généralement robuste) du
test de Student sur deux échantillons indépendants.
On met les deux échantillons ensemble pour former un seul échantillon global de
taille n + p : (U1 , . . . , Un , V1 , . . . , Vp ). On classe ensuite les variables {Ui , Vj } par leur
rang global dans cet échantillon global : cela donne un vecteur de rangs que l’on
note RU,V . On note R1 , . . . , Rn les rangs associés aux variables Ui et S1 , . . . , Sp les
rangs associés aux variables Vj .
R1 = 3, R2 = 5, R3 = 2, S1 = 1, S2 = 4
On pose
Σ1 = R1 + R2 + . . . + Rn , Σ2 = S1 + S2 + . . . + Sp
Principe : pour simplifier, prenons d’abord le cas simple où les deux échantillons sont
de même taille. Alors, sous H0 , on s’attend à ce que Σ1 et Σ2 soit à peu près égaux.
Pour fixer les idées, imaginons que l’alternative corresponde au fait que les Ui ont
tendance à prendre des valeurs supérieures aux Vi . Alors, sous H1 , les rangs Ri des
Ui dans l’échantillon global seront dans l’ensemble supérieurs aux rangs Sj des Vj
dans l’échantillon global. Donc sous H1 , Σ1 sera "grand" (c’est-à-dire "anormalement
grand" par rapport à ce qui se passe sous H0 ).
Maintenant, même si les échantillons ne sont pas de même taille, sous H1 , Σ1
aura tendance à être anormalement grand par rapport à ce qui se passe sous H0 .
Plus généralement, si on pense que U et V n’ont pas la même loi et que l’une des
deux variables a tendance à prendre des valeurs supérieures à l’autre mais on n’a
pas d’intuition sur laquelle des deux, alors s’attend à ce que Σ1 soit "anormalement
grand" ou "anormalement petit" (toujours par rapport à ce qui se passe sous H0 ).
Maintenant la question est : qu’est-ce qu’une valeur "normale" sous H0 ? Les
résultats suivants répondent à cette question.
Proposition 3.23. On a
n(n + 1) n(n + 1)
≤ Σ1 ≤ np +
2 2
p(p + 1) p(p + 1)
≤ Σ2 ≤ np +
2 2
77
Sous H0 : F = G, et sous les conditions 1, 2 et 3 ci-dessus, on a, pour tout i et tout
j,
n+p+1
E(Ri ) = E(Sj ) =
2
(n + p)2 − 1
Var(Ri ) = Var(Sj ) =
12
n(n + p + 1) p(n + p + 1)
E(Σ1 ) = , E(Σ2 ) =
2 2
np(n + p + 1)
Var(Σ1 ) = Var(Σ2 ) =
12
Démonstration.
n(n + 1)
Σ1 ≥ 1 + 2 + . . . + n =
2
et
p(p + 1)
Σ2 ≥ 1 + 2 + . . . + p =
2
Pn+p (n+p)(n+p+1)
Comme Σ1 + Σ2 = i=1 i= 2
, on a
et
1 XN N + 1 2
Var(Ri ) = Var(Sj ) = i2 −
N i=1 2
(N + 1)(2N + 1) N + 1 2
= −
6 2
N2 − 1
=
12
Donc on a
n(n + p + 1)
EΣ1 = E(R1 + . . . Rn ) = nER1 =
2
78
et
p(n + p + 1)
EΣ2 = E(S1 + . . . Sp ) = pES1 =
2
Il reste le calcul des variances de Σ1 et Σ2 . Attention les variables Ri et Sj sont
de même loi mais pas indépendantes ! On a
n
X n X
X
VarΣ1 = Var(R1 + . . . + Rn ) = Var(Ri ) + Cov(Ri , Rj )
i=1 i=1 j6=i
On a déjà calculé les variances, il faut donc calculer maintenant les covariances. Soit
donc i 6= j,
Cov(Ri , Rj ) = E (Ri − ERi )(Rj − ERj )
X N +1 N +1
= (k − )(l − )P(Ri = k, Rj = l)
1≤k,l≤N,k6=l 2 2
Or, (Ri , Rj ) a la même loi que (R1 , R2 ) et, d’après le théorème 7, on a, pour k 6= l,
1
P(R1 = k, R2 = l) =
N (N − 1)
Donc
1 X N +1 N +1
Cov(Ri , Rj ) = Cov(R1 , R2 ) = (k − )(l − )
N (N − 1) k6=l 2 2
Or on a
X X XN
N +1 N +1 N +1 N +1 N +1 2
(k − )(l − )= (k − )(l − )− (k − )
k6=l 2 2 1≤k,l≤N 2 2 k=1 2
X
N 2 N
X
N +1 N +1 2
= (k − ) − (k − )
k=1 2 k=1 2
De plus
N
X XN
N +1 N (N + 1)
(k − )= k− =0
k=1 2 k=1 2
Donc
XN
1 N +1 2 1
Cov(Ri , Rj ) = − (k − ) =− VarR1
N (N − 1) k=1 2 N −1
79
Et finalement
n
X n X
X
Var(Σ1 ) = Var(Ri ) + Cov(Ri , Rj )
i=1 i=1 j6=i
Σ2 ∼ Σ02 (3.1)
80
où Σ02 = S10 + . . . + Sp0 . Or, pour tout j ∈ [p], Sj0 = N + 1 − Sj . Donc
Σ2 ∼ (N + 1)p − Σ2
MV ∼ np − MV = np − (np − MU ) = MU .
np np(n + p + 1)
E(MU ) = Var(MU ) = .
2 12
81
Démonstration. On admet la convergence en loi.
MU = Σ1 − n(n+1) 2
= R1 + . . . + Rn − n(n+1) 2
est une fonction du vecteur
(R1 , . . . , Rn ). On connait la loi de ce vecteur sous H0 , cette loi est donnée par
le théorème 7 : pour toute suite d’entiers (r1 , . . . , rn ) à valeur dans [N ], on a
1
PF =G (R1 , . . . , Rn ) = (r1 , . . . , rn ) =
N (N − 1) . . . (N − n + 1)
Avec R
C’est exactement la même formulation que le test de Wilcoxon des rangs signés,
sauf qu’on met paired=F. C’est en fait False par défaut.
Dans l’exemple lié aux salaires, en supposant cette fois que les échantillons de
salaires d’hommes et de femmes sont i.i.d. et indépendants entre eux, on peut utiliser
wilcox.test(data=salaires,hommes~femmes,alternative="greater")
Remarque 3.28. En plus d’être adaptés à un plus grand nombre de lois, les tests basés sur les
rangs sont plus robustes à la présence d’observations aberrantes, ou "outliers", dans l’échantillon
(penser à la différence médiane/moyenne) .
82
Remarque 3.29. Certains auteurs préconisent, avant l’utilisation éventuelle de Mann-Whitney,
de tester si les deux échantillons ont le même paramètre d’échelle (même variance par exemple). En
effet, si on considère Mann-Whitney comme un test d’égalité en loi supposé détecter une différence
de position, alors il ne semble pas judicieux d’utiliser Mann-Whitney si les échelles diffèrent (ni
d’ailleurs si la forme général de l’histogramme est très différente). Si on fait ce test dans cette
optique-là, alors il parait judicieux de vérifier cette condition sur un graphique par exemple (de toute
façon il faut toujours représenter les données avant toute chose). Il existe aussi des tests d’échelle
(par exemple le test de Levene, qui a des propriétés de robustesse). Citation de Zimmerman (2004) :
"for a wide variety of non-normal distributions, especially skewed distributions, the Type I error
probabilities of both the t test and the Wilcoxon-Mann-Whitney test are substantially inflated by
heterogeneous variances, even when sample sizes are equal."
Cependant, certains praticiens utilisent le test de Mann-Whitney comme un test pour savoir
en gros si l’une des deux populations (U ou V ) a tendance à prendre des valeurs plus grandes que
l’autre. Il n’est alors pas vu comme un test d’égalité en loi. A ce moment-là, on n’a pas besoin
de vérifier si les lois semblent les mêmes à un paramètre de position près (et donc pas besoin de
vérifier que l’échelle est la même).
Remarque 3.30. Une question naturelle : quel type de test (paramétrique/ non
paramétrique) choisir ?
Souvent, si le modèle paramétrique est correct, les tests paramétriques sont plus
puissants que les tests non paramétriques. Cependant, ils sont aussi plus contrai-
gnants, car il faut vérifier les conditions d’application qui sont plus nombreuses dans
ce cas. On choisira généralement un test non paramétrique lorsque
— les conditions d’application du test paramétrique ne sont pas vérifiées
— ou il est impossible de vérifier ces conditions.
"On préconise aussi parfois l’utilisation de tests non paramétriques dans le cas
de petits échantillons, mais le fait d’avoir de petits échantillons ne justifie pas à lui
seul l’utilisation de tests non paramétriques : si les échantillons sont petits, mais
que ce type de données a été suffisamment étudié pour que l’on puisse supposer la
normalité de la distribution, pas de problème pour utiliser des tests paramétriques.
Ce type de conseils est en général donné par prudence, parce que le petit nombre
de données ne permet pas de vérifier, à partir de l’échantillon, la normalité de la
distribution. Dans le doute, on peut donc choisir un test non paramétrique. Les tests
non paramétriques sont certes un peu moins puissants que les tests paramétriques,
mais leur efficacité relative reste bonne" (citation de C. Chabanet, cf biblio).
Remarque 3.31. De la même manière que le test de Student pour comparer les
moyennes de deux échantillons se généralise à plus de deux échantillons par l’ana-
lyse la variance (cf cours de MLG, cas de la régression sur une variable qualitative),
"l’équivalent" du test de Wilcoxon de la somme des rangs pour plus de deux échan-
tillons existe et s’appelle le test de Kruskal-Wallis.
A nouveau, pour Kruskal-Walllis, les données sont remplacées par leur rang dans l’échantillon
global mais cette fois on calcule les sommes de carrés intra-groupe. L’idée est que sous l’hypothèse
nulle (la loi ne dépend pas du groupe) le problème se réduit à nouveau à un problème combinatoire
(il y a une uniformité sous-jacente).
83
Remarque 3.32. Comparons maintenant les tests de Kolmogorov-Smirnov (noté
KS) et le test de Mann-Whitney (noté MW). Le test KS est sensible à tout chan-
gement dans les deux distributions. Des différences substantielles dans la forme,
l’étendue ou la médiane vont amener à une petite p-valeur. En revanche, le test
MW test est seulement sensible à un changement de position (cf plus loin pour une
illustration).
84
Illustration : On regarde ci-dessous les performances respectives des tests de Kolmogorov-Smirnov et de
Mann-Whitney sur un cas particulier. Plus précisément, on utilise deux échantillons qui ne sont pas de
même loi : l’un est de loi normale, l’autre est un mélange de deux lois gamma dont on représente la densité
ci-dessous. On utilise KS et MW pour tester l’égalité des lois sur ces deux échantillons. On regarde la p-valeur
de chaque test.
f=function(x){0.5*(dgamma(x,shape=2,rate=1)+dgamma(-x,2,1))}#densité de probabilité,
#mélange d'une loi gamma et de sa symétrisée
x=seq(-10,10,by=0.1)
plot(x,f(x),type="l")
0.15
0.10
f(x)
0.05
0.00
−10 −5 0 5 10
x
#simulation de 1000 expériences de test et calcul des p-valeurs
KS=rep(0,1000)
MW=rep(0,1000)
for (i in 1:1000){
z=rnorm(100)#simulation de N(0,1)
# simulation d'un échantillon y de loi de densité f
#simulation d'un échantillon t de loi gamma(2,1)
t=c(rgamma(100,shape=2,rate=1))
#simulation de Rademacher
rad=2*rbinom(100,size=1,0.5)-1
y=rad*t
KS[i]=ks.test(y,z)$p.value #p-valeur du test de kolmogorov Smirnov
MW[i]=wilcox.test(y,z)$p.value# pvaleur du test de Mann-Whitney
}
#moyennes des p-valeur de chaque test sur les 1000 simulations
mean(KS)
## [1] 0.002863565
mean(MW)
## [1] 0.4830098
1
Chapitre 4
Formule de Taylor avec reste intégral : on suppose cette fois que f ∈ C n (I) (n
fois continument dérivable) alors, pour tout couple (x, y) de l’interieur de I,
n−1
X
(k) (y − x)k Z y (y − t)n−1 (n)
f (y) = f (x) + f (t)dt
k=0 k! x (n − 1)!
4.2 Introduction
Dans tout le chapitre, l’objectif sera d’estimer une densité f . Pour cela, on s’ap-
puiera sur un n-échantillon iid X = (X1 , . . . , Xn ) où chacune des variables Xi admet
la densité f (par rapport à la mesure de Lebesgue).
Mesure de la qualité d’un estimateur :
1. Définition d’une distance sur l’espace des fonctions :
86
1
R p
— Distance Lp : d(f, g) = kf − gkp = |f (x) − g(x)|p dx
Cas usuel p = 2 ou p = 1.
— distance L∞ : d(f, g) = kf − gk∞ = supx∈R |f (x) − g(x)|.
— Distance ponctuelle en x0 : d(f, g) = |f (x0 ) − g(x0 )|
2. Définition d’une fonction de perte ω : R → R+ telle que ω est convexe et
ω(0) = 0.
Exemple : ω(x) = x2 .
3. Définition du risque d’un estimateur fˆn :
Exemples usuels :
— d(f, g) = |f (x0 ) − g(x0 )|, ω(x) = x2 :
On cherche à déterminer fˆn tel que R(fˆn , f ) soit minimal. Comme expliqué dans
l’introduction, on ne suppose pas que la fonction de densité f appartient à une
famille paramétrique. On va faire une hypothèse moins précise : f appartient à une
classe fonctionnelle qu’on note F. On peut alors définir un risque, qu’on appelle
risque minimax de fˆn sur la classe F , par
On va donc chercher un estimateur fˆn tel que le risque R(fˆn , F) tende vers zéro le
plus vite possible quand n tend vers l’infini.
Définition 4.1. soit (rn )n une suite et une constante C telles que
∀n R(fˆn , F) ≤ Crn
On dit que la suite d’estimateurs (fˆn )n atteint la vitesse (ou le taux) rn sur la classe
F (pour la distance d et la perte ω. )
87
Nous verrons que la vitesse sera d’autant plus grande que la classe F sera une
classe de régularité élevée.
Exemple de classes de fonctions : C k , la classe de Holder (cf définition ci-dessous),
boule dans un espace de Sobolev ( cf cours d’analyse fonctionnelle).
Définition 4.2. Si β ∈ R on note bβc l’entier naturel qui soit le plus grand entier
strictement inférieur à β.
ex : si β = 3, 5 alors bβc = 3 et si β = 4 alors bβc = 3.
Définition 4.3. Pour tout β > 0 et tout L > 0, on définit la classe de Holder de
régularité β et de rayon L par
88
4.3.1 Un estimateur simple de la densité : l’histogramme
Supposons pour simplifier qu’on soit en dimension 1 et que les variables de
l’échantillon soient à valeurs dans [0, 1] donc f : [0, 1] → R+ .
On se donne un découpage de [0, 1] en un certain nombre de classes ]a1 , a2 ], . . . , ]ap , ap+1 ].
Pour simplifier encore, on suppose que les classes sont de même longueur ai+1 − ai =
ai − ai−1 . Cette longueur est notée h. Estimer f par la méthode de l’histogramme
consiste simplement à estimer f par une fonction constante sur chaque classe, cette
constante étant liée à la proportion de Xi tombant dans cette classe. Plus exactement
on pose, pour t ∈]aj , aj+1 ],
1
fˆn (t) = Card{i : Xi ∈]aj , aj+1 ]}
nh
Pour voir très exactement d’où vient cette formule : on a, si f est égale à une
constante cj constante sur ]aj , aj+1 ],
Z aj+1
F (aj+1 ) − F (aj ) = f (t)dt = cj h
aj
89
Code R et illustration graphique du choix du nombre de
classes.
On va illustrer l’importance de bien choisir le nombre de classes par un exemple faisant intervenir une densité
bimodale. On va pour cela simuler un mélange de deux lois gaussiennes : la densité simulée est
1 1 (x − 2)2 (x − 6)2
f (x) = √ exp(− ) + exp(− )
2 2π 2 2
On devrait donc, si l’approximation par l’histogramme est bien faite, se retrouver avec deux “cloches” qui se
chevauchent un petit peu (écart-type=1) et qui sont centrées en 2 et 6 respectivement.
Simulation d’un échantillon de taille n=500 de loi de densité f :
f=function(x){0.5*dnorm(x,mean=2)+0.5*dnorm(x,mean=6)}
sim=function(n){
X=rnorm(n,2,1)
Y=rnorm(n,6,1)
ber=rbinom(n=n,size=1,prob=0.5)
return(ber*X+(1-ber)*Y)}
Z=sim(500)
On estime la densité par un histogramme (on utilise ici la bibliothèque ggplot2) et on rajoute la vraie densité
f en rouge :
library(ggplot2)
p<-ggplot(data.frame(x=Z),aes(x))+labs(x="",y="")
p1<-p+ geom_histogram(aes(y=..density..),color="black",fill="white")+
stat_function(fun=f,col='red')+
labs(title="nb de classes= 30")
p1
0.20
0.15
0.10
0.05
0.00
0.0 2.5 5.0 7.5
La fonction histogram dans ggplot calcule un histogramme avec 30 classes par défaut (ce qu’il signale
d’ailleurs). Ce n’est donc pas la valeur optimale en général. Essayons avec d’autres valeurs du nombre de
classes (=bins).
53
p1<-p+
geom_histogram(aes(y=..density..),bins=3, color="black",fill="white")+
stat_function(fun=f,col='red',xlim=c(-4,12))+
labs(title="nb de classes = 3")
p2<-p+
geom_histogram(aes(y=..density..),bins=10, color="black",fill="white")+
stat_function(fun=f,col='red',xlim=c(-4,12))+
labs(title="nb de classes = 10")
p3<-ggplot(data.frame(x=Z),aes(x))+
geom_histogram(aes(y=..density..),bins=100, color="black",fill="white")+
stat_function(fun=f,col='red',xlim=c(-4,12))+
labs(title="nb de classes = 100",x="",y="")
0.15
0.15
0.2
0.10
0.10
0.1
0.05 0.05
p1<-p+
geom_histogram(aes(y=..density..),color="black",fill="white")+
stat_function(fun=f,col='red')+
labs(title="nb de classes=30")
p2<-ggplot(data.frame(x=Z),aes(x))+
54
geom_histogram(aes(y=..density..),bins=100,color="black",fill="white")+
stat_function(fun=f,col='red')+
labs(title="nb de classes=100")
grid.arrange(p1,p2,nrow=1)
0.15 0.15
density
0.10 0.10
0.05 0.05
0.00 0.00
0 4 8 −2.5 0.0 2.5 5.0 7.5 10.0
x
On voit donc qu’avec un nombre de classes égal à 100, on a, conrairement à précédemment, un très bon choix.
La taille optimale du nombre de classes est croissante avec n, autrement dit, le pas h optimal décroit avedc n,
ce que l’on va illustrer plus tard avec l’estimateur à noyau de fenêtre h.
Remarquez que l’on fait deux approximation successives : une première approximation quand on approche la
densité par une fonction constante par morceaux, et ensuite une deuxième approximation quand on approche
chaque constante à l’aide des données.
55
4.3.2 Estimateurs à noyaux
Un inconvénient de l’estimateur par histogramme précédent est que la fonction
de densité résultante fˆn n’est pas régulière : il s’agit d’une fonction constante par
morceau, qui a donc des sauts aux extrémités de chaque classe. En général, la densité
à estimer est plus lisse, au moins continue.
L’estimation par noyau a pour but de répondre à cet écueil.
Principe : Si f est continue en x (ce qui va être le cas pour les classes de fonctions
qu’on va considérer) alors
F (x + h) − F (x) F (x + h) − F (x − h)
f (x) = F 0 (x) = lim = lim
h→0 h h→0 2h
L’idée est donc d’utiliser l’approximation suivante, pour h petit,
F (x + h) − F (x − h)
f (x) ≈
2h
Pour estimer la densité f on peut donc passer par un estimateur F̂n de la cdf F .
Voyons ce qui se passe si on choisit comme estimateur la fonction de répartition
P
empirique Fn . (On rappelle que Fn (x) = n1 ni=1 1Xi ≤n ) On choisit un h > 0 petit
pour que l’approximation ci-dessus soit valable, et on pose
n
Fn (x + h) − Fn (x − h) 1X 1
f˜n (x) = = 1X ∈]x−h,x+h]
2h n i=1 2h i
93
Définition 4.7. Etant donné K un noyau et h > 0, on pose
n
1X 1 Xi − x
∀x ∈ R, fˆn (x) = K( )
n i=1 h h
Remarque 4.8. — La plupart des noyaux sont symétriques, positifs et sont dé-
croissants sur R+ comme le noyau Gaussien : plus y est proche de 0, plus
K(y) est grand. Donc, pour un x ∈ R donné, plus une observation Xi est
proche de x, plus K( Xih−x ) est grand. Donc fˆn (x) est d’autant plus grand que
x est proche de beaucoup d’observations Xi (somme de beaucoup de grandes
valeurs K( Xih−x )).
— L’estimateur est somme de fonctions K( Xih−x ) qui sont continues si K est
continu. Donc fˆn est continu si K est continu.
R
— fˆn (x)dx = 1, donc, si K(x) ≥ 0 ∀x ∈ R, alors fˆn est une densité.
— Le paramètre h > 0 est appelé fenêtre (bandwidth). C’est un paramètre de
lissage : plus h est grand, plus l’estimateur est régulier. Comme dans le cas
de l’estimateur à histogramme, le choix de h est délicat, la fenêtre h optimale
devant réaliser un équilibre biais/variance (cf section suivante).
— Dans la pratique, le choix du noyau est peu influent, contrairement au choix
de la fenêtre !
94
Illustration graphique et code R
On va utiliser le même exemple de distribution bimodale que précédemment. L’estimation par noyaux peut
se faire avec différentes méthodes. On peut utiliser la fonction density du package stat. Cette procédure
n’estime que des densités à une seule variable. Pour des fonctions multivariées, on peut utiliser par exemple
la fonction kde du package ks (de 1 à 6 variables).
Par défaut le noyau utilisé est le noyau gaussien, il est possible de changer de noyau avec l’option kernel.
On va en fait utiliser la version de ggplot pour représenter l’estimateur à noyau. La fonction qui permet de
dessiner l’estimateur à noyau est
geom_density
Le paramètre représentant le fenêtre h s’appelle bw (comme bandwidth).
On illustre l’influence du choix de la fenêtre. On tire les mêmes conclusions que pour l’histogramme.
p<-ggplot(data.frame(x=Z),aes(x))+labs(x="",y="")
p1<-p+geom_density(bw=0.1)+stat_function(fun=f,col='red',alpha=0.4)+ggtitle("h=0.1")
p2<-p+geom_density(bw=0.5)+stat_function(fun=f,col='red',alpha=0.4)+ggtitle("h=0.5")
p3<-p+geom_density(bw=0.8)+stat_function(fun=f,col='red',alpha=0.4)+ggtitle("h=0.8")
p4<-p+geom_density(bw=1.2)+stat_function(fun=f,col='red',alpha=0.4)+ggtitle("h=1.2")
grid.arrange(p1,p2,p3,p4,nrow=2,ncol=2)
h=0.1 h=0.5
0.25 0.20
0.20
0.15
0.15
0.10
0.10
0.05
0.05
0.00 0.00
0.0 2.5 5.0 7.5 0.0 2.5 5.0 7.5
h=0.8 h=1.2
0.20 0.20
0.15 0.15
0.10 0.10
0.05 0.05
0.00 0.00
0.0 2.5 5.0 7.5 0.0 2.5 5.0 7.5
Pour finir, on illustre le choix de deux fenêtres calculées à partir des données. L’une est la méthode de Sheather
et Jones (SJ) et l’autre est basée sur la validation croisée, qui sera vue en fin de chapitre (ucv=unbiased
58
cross-validation). Pour d’autres méthodes, consultez la documentation liée à bw.
p<-ggplot(data.frame(x=Z),aes(x))+labs(x="",y="")
p5<-p+geom_density(bw="ucv")+stat_function(fun=f,col='red',alpha=0.4)+ggtitle("ucv")
p6<-p+geom_density(bw="SJ")+stat_function(fun=f,col='red',alpha=0.4)+ggtitle("SJ")
grid.arrange(p5,p6,ncol=2)
ucv SJ
0.20 0.20
0.15 0.15
0.10 0.10
0.05 0.05
0.00 0.00
0.0 2.5 5.0 7.5 0.0 2.5 5.0 7.5
Il existe une version en dimension 2 de cette fonction dans ggplot2 qui s’appelle
geom_density_2d
.
59
4.4 Risque quadratique ponctuel des estimateurs
à noyau sur la classe des espaces de Holder
Dans cette section, on s’intéresse au risque quadratique ponctuel de fˆn , i.e. étant
donné x0 ∈ R
R(fˆn , f ) = E |fˆn (x0 ) − f (x0 )|2
Définition 4.9. Soit ` ∈ N∗ . ROn dit que le noyau K est d’ordre ` si ∀j ∈ {1, . . . , `},
u → uj K(u) est intégrable et uj K(u)du = 0.
Proposition 4.10.R
Si f ∈ Σ(β, L) avec β > 0 et L > 0 et si K est un noyau d’ordre
` = bβc tel que |u|β |K(u)|du < ∞ alors pour tout x0 ∈ R, et pour tout h > 0 le
biais peut être borné comme suit :
hβ L Z
|E[fˆn (x0 )] − f (x0 )| ≤ |u|β |K(u)|du
`!
Démonstration. On a
h1 X
1 n
X i − x0 i
E[fˆn (x0 )] = E K( )
n i=1 h h
h1 X 1 − x0 i
E K( )
h Z h
1 u − x0
= K( )f (u)du
Zh h
= K(v)f (x0 + hv)dv
De plus Z
f (x0 ) = f (x0 ) × 1 = f (x0 ) K(v)dv.
Donc h i Z h i
E fˆn (x0 ) − f (x0 ) = K(v) f (x0 + hv) − f (x0 ) dv
97
pour un certain ξ ∈]0, 1[.Donc
Z h i Z hX
`−1
(hv)k (k) (`) (hv)` i
K(v) f (x0 + hv) − f (x0 ) dv = K(v) f (x0 ) + f (x0 + hvξ) dv
k=1 k! `!
h` Z
= K(v)v ` f (`) (x0 + hvξ)dv
`!
R
Comme K est d’ordre `, on a aussi K(v)v ` f (`) (x0 )dv = 0. Donc on a
Z h i h` Z h i
K(v) f (x0 + hv) − f (x0 ) dv = K(v)v ` f (`) (x0 + hvξ) − f (`) (x0 ) dv
`!
Or, comme f ∈ Σ(β, L), on a |f (`) (x0 + hvξ) − f (`) (x0 )| ≤ L|hv|β−` . Et finalement
Z
h i |h|` Z
K(v) f (x0 + hv) − f (x0 ) dv ≤ |K(v)||v|` L|hv|β−` dv
`!
ce qui signifie que
h
i
E fˆn (x0 ) − f (x0 )
L|h|β Z
≤ |K(v)||v|β dv
`!
Le biais au carré tend donc vers zéro à la vitesse h2β . Plus la fonction f est
régulière, plus le biais tend vite vers zéro quand h tend vers zéro (à condition bien
sûr que l’ordre du noyau soit suffisamment grand).
Proposition 4.11. Si f est bornée et si K est de carré intégrable alors
kf k∞ kKk22
Var(fˆn (x0 )) ≤
nh
En particulier, si f ∈ Σ(β, L) alors
M (β, L)kKk22
Var(fˆn (x0 )) ≤
nh
Démonstration.
1 Xn
X i − x0
Var(fˆn (x0 ) = Var K( )
nh i=1 h
n
X 1 X i − x0
= Var K( )
i=1 nh h
n
X 1 X i − x0
= 2 2
Var K( )
i=1 n h h
1 X 1 − x0
= Var K( )
nh2 h
1 2 X 1 − x0
≤ E K ( )
nh2 Z h
1 2 u − x0
= K f (u)du
nh2Z h
1
= K 2 (v)f (x0 + vh)dv
nh
98
Et enfin, on utilise la proposition 10 : il existe une constante positive M (β, L) tel
que kf k∞ ≤ M (β, L). Ceci implique que
Z
ˆ 1
Var(fn (x0 ) ≤ M (β, L) K 2 (v)dv
nh
Pour que la variance tende vers zéro, il faut que nh tende vers l’infini. En parti-
culier, à n fixé, la variance est une fonction décroissante de h au contraire du biais
qui est une fonction croissante de h. Il y a donc une valeur optimale de h qui doit
réaliser l’équilibre entre le biais au carré et la variance. On peut à présent donner
un contrôle du risque quadratique.
Théorème R4.12. Soit β > 0 et L > 0 et K un noyau de carré intégrable et d’ordre
bβc tel que |u|β |K(u)|du < ∞. Alors, en choisissant une fenêtre de la forme h =
1
cn− 2β+1 avec une constante c > 0, on obtient
2β
∀x0 ∈ R, R(fˆn (x0 ), Σd (β, L) := sup E[|fˆn (x0 ) − f (x0 )|2 ] ≤ Cn− 2β+1
f ∈Σd (β,L)
99
1
Autrement dit, pour une fenêtre h de l’ordre de n− 2β+1 , le biais au carré et la variance
1
sont de même ordre. Plus exactement, si on choisit la fenêtre h∗ = cn− 2β+1 , avec c
une constante positive, on a
1
Biais au carré ≈ h2β
∗ ≈ variance ≈
nh∗
De plus on a alors
−2β
h2β
∗ ≈ n
2β+1
Autrement dit, il existe une certaine constante C telle que, pour cette fenêtre h∗ , on
a
−2β
ˆ
R fn (x0 ), Σd (β, L) ≤ Cn 2β+1
Cette fenêtre est donc optimale à une constante près (si on change c, on change C
−2β
mais ça ne change pas le taux qui est n 2β+1 ).
100
R1
— −1 φm (u)φk (u)du = 1m=k
— φm est un polynôme de degré m.
— φ2m est pair et φ2m+1 est impair ∀m ≥ 0.
P`
Proposition 4.14. Soit K` : u → m=0 φm (0)φm (u)1|u|≤1 . Alors K` est un noyau
d’ordre `.
Donc
Z Z 1 X
j
uj K(u)du = aq φq (u)K(u)du
−1 q=0
j
X Z 1 X̀
= aq φq (u) φm (0)φm (u)du
q=0 −1 m=0
j X̀
X Z 1
= aq φm (0) φq (u)φm (u)du
q=0 m=0 −1
j
X
= aq φq (0)
q=0
0 si j ≥ 1
=
1 si j = 0
Remarque 4.15. Comme φ2k+1 est impaire, on a φ2k+1 (0) = 0 et donc K2k =
K2k+1 . Et donc l’ordre maximal de K` est impair.
Or la fonction f étant inconnue, ce risque n’est pas calculable à partir des don-
nées. On cherche donc à estimer ce risque en utilisant uniquement les données.
Remarquons tout de suite que minimiser en h la quantité R(fˆn,h , f ) est équivalent à
101
minimiser en h la quantité R(fˆn,h , f ) − kf k22 . On va en fait remplacer la minimisa-
tion de la quantité inconnue R(fˆn,h , f ) − kf k22 par la minimisation d’un estimateur
R̂(h) de cette quantité. Plus précisément on va chercher un estimateur sans biais de
R(fˆn,h , f ) − kf k22 .
Pour simplifier on suppose dans le théorème R
suivant que K est positif (on aurait
pu aussi supposer que f et K sont tels que |K( u−v h
)|f (u)f (v)dudv est finie). De
cette manière toutes les quantités que l’on manipulera seront positives (car K et f
sont positives) et on pourra appliquer Fubini. On suppose aussi que R(fˆn,h , f ) < ∞
et f ∈ L2 .
2 X X n
1 X i − Xj
R̂(h) = kfˆn,h k22 − K
n(n − 1) i=1 j=1,j6=i h h
Or
Z
R(fˆn,h , f ) − kf k22 = E kfˆn,h k22 − 2 fˆn,h (x)f (x)dx
Z
= Ekfˆn,h k22 − 2 Efˆn,h (x)f (x)dx
102
On définit alors
ĥ = arg min R̂(h)
h∈H
si ce minimum est atteint. On cherche une fenêtre parmi une grille finie de valeurs,
grille qu’on a notée H dans la formule ci-dessus.
L’estimateur fˆn,ĥ a de bonnes propriétés pratiques et des propriétés de consis-
tance.
La validation croisée est une méthode très générale dont on reparlera plus en
détail dans le prochain chapitre. L’idée d’utiliser un estimateur sans biais du risque
est aussi une idée assez générale (cf critère Cp).
103
Chapitre 5
5.1 Introduction
Dans ce chapitre, on cherche à expliquer les valeurs que peut prendre une variable
Y à partir des valeurs que peut prendre une variable X.
Exemples :
— Y est le taux d’insuline dans le sang, qu’on explique (ou prédit) à l’aide de
X= (IMC, pression du sang, concentration de molécules).
— Y est le niveau de diplôme obtenu , qu’on explique à l’aide de X = (âge, sexe,
revenu des parents, métier des parents).
On suppose que la variable Y est intégrable E|Y | < ∞ et on note r la fonction
de régression de Y sur X :
r(x) = E(Y |X = x)
L’objectif est d’estimer la fonction r pour expliquer et prédire Y à partir de X. Pour
cela on dispose des réalisations de n couples de variables (X1 , Y1 ), . . . , (Xn , Yn ). On
va supposer que les (Xi , Yi ) sont indépendants.
Vocabulaire
— les Yi sont les variables à expliquer ou les variables réponses ou variables de
sortie.
— les Xi constituent le design, les variables explicatives, les covariables, ou va-
riables d’entrée.
Modélisation
Le design pourra être aléatoire ou déterministe. Dans ce dernier cas, on notera
plutôt xi à la place de Xi .
Le fait que r(x) = E(Y |X = x) se réécrit
104
Les i sont appelées erreurs et jouent le rôle de bruit. Dans la suite, on va faire
une hypothèse très forte :
r(x) = β0 + β1 x1 + . . . , βp xp
pour x ∈ Rp .
105
●
●
●
●
●● ● ●
● ● ●●● ●
● ●● ●●
●●● ●●●●●●●●
25
●●● ●
●●●●●●●●●●●●●●
● ● ●●● ●● ●●
● ● ●●●●●●●●●●●●●●●
●●●●● ●●
●●●●●●●●●●●●●●● ●
● ●●●●●●● ● ● ●
● ●●●●●●●●●●●●●●●● ●
●●● ●●●●●●●●●●●●●●●
●●●●●●●●●●●●●●●●●● ●
● ●● ●●● ●●●●●●●●●●● ●
●●●●●●●●●●●●●●●●●●
●● ●●●●●●●●●●●●●●●● ●
●●●●●●●●●●●●●●●●●●●●
●●●●●●●●●● ●●● ● ● ●
●●●●●●●●●●● ●●● ●●●
● ●●●●●●●●●● ●● ● ●
● ●● ●●●●●●●●●●●●● ●
●●● ●●●●●●●●●●●●●● ●
●● ●●●●●●●●●●●●●● ●
20
● ●●●●●●●●●●● ●●●●
●●●●●●●●●●●●●●●●●● ●
ht
●●●●●●●●●●● ●●●
● ●●●●●●●●●●● ●
● ●●● ●●●●●●● ●● ●
●●●●●● ● ●
●●● ●●●●●●● ●
● ●●● ●●● ●●● ●●●
●●● ●●●●● ●
● ●●●●●●● ●●●
●● ●●● ●● ●
● ● ●●●● ●●
● ●●● ●● ●●
●●●●●●●● ●
● ●●● ●
● ● ● ●● ●●
●●● ●●●●●● ●
● ●●● ●
● ● ●
●●●● ●
15
● ●● ●
●●
●● ● ●
● ●
●●●
● ●
●● ●
●
●
● ●
●
●
30 40 50 60 70
circ
Figure 5.1 – Représentation hauteur versus circonférence pour les 1429 eucalyptus
mesurés
Si les données se trouvent dans un data.frame appelé euca et si les noms des
variables sont ht et circ alors on peut utiliser
reg=lm(ht~circ,data=euca)
On peut ensuite représenter le nuage de points avec la droite de régression, ainsi que
l’intervalle de confiance sur un ensemble de valeurs de prévisions (à 95%) .
> plot(ht~circ,data=euca)
> circ=euca[,’circ’]
> grille<-seq(min(circ),max(circ),length=100)
> grilledataframe<-data.frame(circ=grille)
> ICpred<-predict(reg,new=grilledataframe,interval="pred",level=0.95)
> matlines(grille,ICpred,lty=c(1,2,2),col=c(’red’,’blue’,’blue’))
Nous constatons que les observations sont globalement bien ajustées par le mo-
dèle, sauf peut-être pour les faibles valeurs de circonférences, qui semblent en ma-
jorité situées en dessous de la droite. Ceci suggère d’utiliser plutôt le modèle de
régression suivant √
ht = a1 + a2 circ + a3 circ +
On peut donc utiliser un modèle linéaire avec une transformation de
la variable d’origine. On peut d’ailleurs vérifier qu’en introduisant la variable
sqrt(circ), on a bien un meilleur modèle :
> reg1=lm(ht~circ,data=euca)
> reg2=lm(ht~circ+I(sqrt(circ)),data=euca)
> anova(reg1,reg2)
106
●
●
●
●
●● ● ●
● ● ●●● ●
● ●● ●●
●●● ●●●●●●●●
25
●●● ●
●●●●●●●●●●●●●●
● ● ●●● ●● ●●
● ● ●●●●●●●●●●●●●●●
●●●●● ●●
●●●●●●●●●●●●●●● ●
● ●●●●●●● ● ● ●
● ●●●●●●●●●●●●●●●● ●
●●● ●●●●●●●●●●●●●●●
●●●●●●●●●●●●●●●●●● ●
● ●● ●●● ●●●●●●●●●●● ●
●●●●●●●●●●●●●●●●●●
●● ●●●●●●●●●●●●●●●● ●
●●●●●●●●●●●●●●●●●●●●
●●●●●●●●●● ●●● ● ● ●
●●●●●●●●●●● ●●● ●●●
● ●●●●●●●●●● ●● ● ●
● ●● ●●●●●●●●●●●●● ●
●●● ●●●●●●●●●●●●●● ●
●● ●●●●●●●●●●●●●● ●
20
● ●●●●●●●●●●● ●●●●
●●●●●●●●●●●●●●●●●● ●
ht
●●●●●●●●●●● ●●●
● ●●●●●●●●●●● ●
● ●●● ●●●●●●● ●● ●
●●●●●● ● ●
●●● ●●●●●●● ●
● ●●● ●●● ●●● ●●●
●●● ●●●●● ●
● ●●●●●●● ●●●
●● ●●● ●● ●
● ● ●●●● ●●
● ●●● ●● ●●
●●●●●●●● ●
● ●●● ●
● ● ● ●● ●●
●●● ●●●●●● ●
● ●●● ●
● ● ●
●●●● ●
15
● ●● ●
●●
●● ● ●
● ●
●●●
● ●
●● ●
●
●
● ●
●
●
30 40 50 60 70
circ
Model 1: ht ~ circ
Model 2: ht ~ circ + I(sqrt(circ))
Res.Df RSS Df Sum of Sq F Pr(>F)
1 1427 2052.1
2 1426 1840.7 1 211.43 163.8 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
107
\color{red}
lm(y~x+x^2))$coefficients\color{black}
lm(y~poly(x,2))$coefficients
108
On peut utiliser les estimateurs à noyau du chapitre précédent :
n
1 X Xi − x
fˆn,X (x) = K( )
nh i=1 h
n
ˆ 1 X Xi − x Yi − y
fn (x, y) = 2
K( )K( )
nh i=1 h h
Proposition 5.1. Si K est un noyau d’ordre 1 alors ∀x ∈ R
Pn X −x
Pi=1 Yi K( ih ) Pn
n X −x i
si i=1 K( Xih−x ) 6= 0
r̂n (x) = i=1
K( h
)
0 sinon
P
Démonstration. fˆn,X (x) = 0 est équivalent à ni=1 K( Xih−x ) = 0.
P
Supposons donc que ni=1 K( Xih−x ) 6= 0. Alors
Z
y fˆn (x, y)
r̂n (x) = dy
fˆn,X (x)
Z n
1 1 X Xi − x Yi − y
= y K( )K( )dy
fˆn,X (x) nh2 i=1 h h
nh 1 X n
Xi − x Z Yi − y
= Pn Xi −x 2
K( ) yK( )dy
i=1 K( h ) nh i=1 h h
1 n
X Xi − x 1 Z Yi − y
= Pn Xi −x K( ) yK( )dy
i=1 K( h ) i=1 h h h
n
X
1 Xi − x
= Pn Xi −x K( )Yi
i=1 K( h ) i=1 h
109
La fenêtre optimale équilibrant biais (au carré) et variance se trouve entre ces
deux extrêmes.
Remarque 5.5. Il se peut que la densité fX soit connue. Dans ce cas, il est préfé-
rable d’utiliser R
yfˆn (x,y) dy si f (x) 6= 0
fX (x) X
r̃n (x) =
0 si f (x) = 0
X
où Z Z
2
y M (x, y)dy < ∞ et y 2 f (x, y)dy < ∞
110
Alors, si |h| ≤ , il existe une constante C(x) (dépendant de x) telle que
1
E[(r̃n (x) − r(x))2 ] ≤ C(x)(h2 + )
nh
De plus si on choisit une fenêtre h telle que h n−1/3 (le signe signifie “de l’ordre
de”), il existe une constante C 0 (x) telle que
— Biais
On va prouver dans le calcul de la variance que,
h sous les hypothèses
i de l’énoncé,
Var[Y1 K( X1h−x )] < ∞. Ceci implique que E |Y1 K( X1h−x )| < ∞. On va pou-
voir utiliser le théorème du transfert pour calculer cette intégrale ainsi que le
théorème de Fubini si besoin.
On a
Xn
1 Xi − x
E[r̃n (x)] = E Yi K( )
nhfX (x) i=1 h
1 X1 − x
= E[Y1 K( )]
hfX (x) h
Z Z
1 t−x
= yK( )f (t, y)dtdy
hfX (x) h
1 Z Z
= yK(v)f (x + vh, y)dvdy
fX (x)
De plus
r(x) = E[Y |X = x]
Z
= yfY (y|X = x)dy
Z
f (x, y)
= y dy
fX (x)
Donc Z
r(x)fX (x) = yf (x, y)dy
Donc on a aussi
Z
r(x + vh)fX (x + vh) = yf (x + vh, y)dy
111
Donc
On a utilisé le fait que K est à support dans [−1, 1] dans la dernière égalité.
On applique l’inégalité des accroissements finis à r et fX car elles sont conti-
nument dérivables au voisinage de x. Il existe une constante C(x) telle que,
pour tout |u| ≤ ,
|r(x + u) − r(x)| ≤ C(x)u
|fX (x + u) − fX (x)| ≤ C(x)u
On peut donc appliquer ces inégalités avec u = vh pour |v| ≤ 1 et |h| ≤ , ce
qui donne
De plus r étant continue sur [x − , x + ], il existe une constante c(x) telle que
|r(x + hv)| ≤ c(x) pour tout |h| ≤ et tout |v| ≤ 1. Donc on a
112
Xn
1 Xi − x
Var(r̃n (x)) = Var Yi K( )
nhfX (x) i=1 h
1 X1 − x
= nVar Y1 K( )
nhfX (x) h
1 X1 − x
=n 2 2 2 Var Y1 K( )
n h fX (x) h
1 2 2 X1 − x
≤ E Y 1 K ( )
nh2 fX2 (x) h
Z
1 t−x
= 2 2
y2K 2( )f (t, y)dtdy
nh fX (x) h
Z
1
= 2
y 2 K 2 (v)f (x + vh, y)dvdy
nhfX (x)
Comme |h| ≤ , on a |hv| ≤ pour tout v ∈ [−1, 1]. Donc , d’après la troisième
hypothèse de l’énoncé,
|f (x + hv, y) − f (x, y)| ≤ M (x, y)
et donc
f (x + hv, y) ≤ f (x, y) + M (x, y) (5.1)
Ainsi
Z Z
1 2 2
Var(r̃n (x)) ≤ y K (v)M (x, y)dvdy + y 2 K 2 (v)f (x, y)dvdy
nhfX2 (x)
R 2
K (v) Z 2 Z
2
= y M (x, y)dy + y f (x, y)dvdy
nhfX2 (x)
113
L’estimateur de Nadaraya-Watson est un cas particulier des estimateurs par
polynomes locaux.
Démonstration.
r̂n (x) = arg min τ (θ)
θ∈R
où n
X X − x
i
τ (θ) = K (Yi − θ)2
i=1 h
τ est un polynôme du second degré en θ. Recherche d’un point critique :
n
X X − x Xn X − x
0 i i
τ (θ) = 0 ⇔ K Yi = θ K
i=1 h i=1 h
Pn
Xi −x
i=1
h
K Yi
⇔θ= Pn
Xi −x
i=1 K h
Pn
Xi −x
C’est un minimum car τ 00 ≡ 2 i=1 K h
≥ 0.
L’estimateur par polynômes locaux est une généralisation de l’estimateur de
Nadaraya-Watson associée à sa caractérisation par la proposition précédente. Il faut
garder à l’esprit ici que l’idée est de regarder les choses localement, et donc que x
est fixé. On aura donc calculé pour ce x fixé un estimateur de r(x) mais si on veut
r̂(y) il faut faire un autre calcul.
L’idée associée à l’estimateur par polynômes locaux est de reprendre le problème
de minimisation de la proposition précédente mais au lieu d’utiliser une constante
θ, on utilise un polynôme.
Plus précisément, si r est régulière alors, autour de x, r est proche du polynôme
associé à son développement de Taylor-Lagrange en x : pour u proche de x on a
avec
X̀ r (k) (x)
P`,x (u) = (u − x)k
k=0 k!
114
Evidemment P`,x est tout aussi inconnu que r(x) (ses coefficients dépendent de la
quantité que l’on cherche à estimer r(x) mais aussi des dérivées r0 (x), . . . , r(`) (x)).
On va en fait essayer d’estimer ce polynôme P`,x . Si on écrit
Comme on n’a pas accès à r(Xi ) mais à sa donnée bruitée Yi , on cherche en fait P̂
tel que
P̂ (Xi ) est proche de Yi pour les Xi proches de x.
Autrement dit
(P̂ (Xi ) − Yi )2 petit pour les Xi proches de x.
Des poids K( Xih−x ) sont ajoutés pour prendre en compte cette notion de proxi-
mité. On pose alors
On pose θ̂ = (θ̂0 , θ̂1 , . . . , θ̂` ). L’estimateur par polynôme local d’ordre ` est alors défini
par
r̂n` (x) = θ̂0
115
On peut aussi écrire r̂(x) = ω(x)T Y où Y est le vecteur (Y1 , . . . , Yn )T et ω(x) =
(ω1 (x), . . . , ωn (x))T .
On a vu que l’estimateur de Nadaraya-Watson est linéaire.
Attention : ne pas confondre le fait que l’estimateur soit linéaire, ce qui sous
entend linéaire en Y , et le fait que la fonction de régression soit linéaire, ce qui
signifie que r(x) est linéaire en x (et on cherche alors un estimateur linéaire en
x). L’estimateur associé aux MCO r̂(x) = β̂ T x est linéaire en x et c’est également
un estimateur linéaire : r̂(x) = xT β̂ = xT (X T X)−1 X T Y = ω(x)T Y où ω(x) =
[xT (X T X)−1 X T ]T est un vecteur qui ne dépend pas de Y .
Introduisons, pour la proposition suivante, quelques notations : pour tout i =
1, . . . , n et tout u ∈ R,
1
u
Xi − x
Zi = , V` (u) =
..
h .
u`
`!
Et on pose
n
X
Bn,x = K(Zi )V` (Zi )V` (Zi )T .
i=1
Proposition 5.11. Si la matrice Bn,x est définie positive alors l’estimateur par
polynômes locaux r̂n` (x) est un estimateur linéaire.
Démonstration. On a
r̂n,` (x) = θ̂0 (x) = eT1 θ̂(x)
avec
1
0
e1 = .
.
.
0
θ̂(x) = arg min
`+1
τ (θ)
θ∈R
où n
X
τ (θ) = K(Zi )(Yi − θT V` (Zi ))2
i=1
On a
n
X h i
τ (θ) = K(Zi ) Yi2 + (θT V` (Zi ))2 − 2Yi θT V` (Zi )
i=1
n
X n
X n
X
= K(Zi )Yi2 + K(Zi )θT V` (Zi )V` (Zi )T θ − 2θT K(Zi )Yi V` (Zi )
i=1 i=1 i=1
= a + θT Bn,x θ − 2θ b T
P Pn
avec a = ni=1 K(Zi )Yi2 et b = i=1 K(Zi )Yi V` (Zi )
Rappels :
116
— Si f (x) = xT a alors ∇f (x) = a et Hf (x) = 0 (Hf est la hessienne de f ).
— Si f (x) = xT Ax alors ∇f (x) = (A + AT )x et Hf (x) = A + AT
— Si A est symétrique et f (x) = xT Ax alors ∇f (x) = 2Ax et Hf (x) = 2A
Recherche de point critique :
Donc
∇τ (θ) = 0 ⇔ Bn,x θ = b
Si Bn,x est définie positive, elle est inversible et donc il y a un seul point critique
donné par
−1
θ̂ = Bn,x b
Ce point critique correspond bien à un minimum global car la fonction est convexe.
En effet
Hτ (θ) = 2Bn,x > 0
On a donc
−1
r̂n,` (x) = eT1 Bn,x b
hX
n i
−1
= eT1 Bn,x K(Zi )Yi V` (Zi )
i=1
n
X
= ωi (x)Yi
i=1
avec
−1
ωi (x) = K(Zi )eT1 Bn,x V` (Zi )
ωi (x) ne dépend que de x, K, `, h, et des Xi et pas des Yi . Donc r̂n,` est bien un
estimateur linéaire.
Remarque 5.12. On a
n
X
ωi (x) = 1
i=1
117
On note maintenant rh l’estimateur utilisant la fenêtre h. Si on enlève une partie
de l’échantillon (Xi , Yi )i∈I avec I une partie de {1, . . . , n} on notera r̂h−I l’estimateur
calculé à partir de l’échantillon auquel on a ôté (Xi , Yi )i∈I .
Remarquez que la fonction de régression r est telle que
h i
r = arg min E (Y − f (X))2 .
f ∈L2 (PX )
On ne peut pas minimiser ce risque puisque r est inconnu. Une première idée est
de remplacer r(Xi ) par son observation bruitée Yi et d’oublier l’espérance, c’est-à-
dire de minimiser n
1X
R̂n (h) = (r̂h (Xi ) − Yi )2
n i=1
NB : cette quantité est connue sous le nom de "erreur d’apprentissage" (training
error).
C’est en général une très mauvaise idée d’utiliser ce risque comme substitut du
vrai risque pour la sélection de modèle ! En effet les mêmes données sont utilisées à
la fois pour estimer r et estimer le risque. Il y a un manque d’indépendance.
Prenons l’exemple de l’EMC non paramétrique. Imaginons qu’on cherche à ajus-
ter un polynôme. On se pose donc la question du degré M . Pour chaque M on
calcule β̂ M l’EMC associé au design X = (Xij )1≤j≤M,1≤i≤n avec Xij = xj−1 i . Si M
est assez grand et si les points du design sont distincts alors le risque empirique est
égal à 0. On a obtenu un polynôme qui passe par tous les points (Xi , Yi ) ("on recopie
les données"). Mais la variance de cet estimateur risque fort d’être trop grande.
L’erreur d’apprentissage est trop optimiste. On aura en général E[R̂n (h)] < R(h).
Utiliser cette erreur pousse au sur-ajustement (overfitting) : l’estimateur associé sera
trop adapté aux données particulières qu’on a et ne se généralisera pas bien à de
nouvelles données.
iid
Remarque 5.14. — Si Y1 , . . . , Yn ∼ Y alors pour estimer E(Y ) on utilise sou-
P
vent son équivalent empirique n1 ni=1 Yi .
— Si g est une fonction fixe (i.e. ne dépendant pas des données) alors Yi −
iid P
g(Xi ) ∼ Y − g(X). Et il est alors naturel d’utiliser n1 ni=1 (Yi − g(Xi ))2 pour
estimer E(Y − g(X))2 . En effet si g est fixe,
n
1X
E (g(Xi ) − Yi )2 = kg − rk2L2 (PX ) + σ 2 ,
n i=1
et n
1X 1
Var (g(Xi ) − Yi )2 = Var(g(X) − Y )2 .
n i=1 n
Si on se donne un ensemble de fonctions déterministes (gh )h∈H dépendant
d’un paramètre h (on entend par "déterministe" le fait que gh ne dépend pas
de l’échantillon), alors minimiser le risque empirique semble un bon substitut à
la minimisation du risque quadratique kgh − rk2L2 (PX ) pour choisir le paramètre
h.
118
5.5.2 Validation croisée
La technique de validation croisée est très générale et s’applique à de nombreuses
procédures d’estimation. Ici on va l’appliquer pour le choix de la fenêtre h de l’es-
timateur par polynômes locaux, mais elle aurait pu être utilisée pour le choix d’un
autre paramètre d’ajustement (le degré du polynôme si on ajuste un polynôme par
les moindres carrés par exemple).
On se donne une grille de valeurs H de fenêtres, parmi lesquelles on veut choisir
une fenêtre optimale ĥ en se basant sur les données uniquement.
Le principe général est de diviser l’échantillon en un ensemble d’apprentissage
(training set) et un ensemble de validation (validation set). On fabrique des estima-
teurs à partir de l’ensemble d’apprentissage et ensuite l’ensemble de validation est
utilisé pour estimer leur risque de prédiction. Les schémas les plus populaires sont
les suivants :
— Hold-out CV : on divise l’échantillon en deux parties I1 et I2 (I1 et I2 sont donc
deux ensembles disjoints de {1, . . . , n}). On calcule les estimateurs (r̂hI1 )h∈H à
partir de (Xi , Yi )i∈I1 . Puis on calcule les estimateurs des risques associés
1 X
R̂(h) = (Yi − r̂hI1 (Xi ))2
n2 i∈I2
où on a noté n2 = Card(I2 ).
— V -fold CV : les données sont divisées en V ensembles disjoints I1 , . . . , IV .
Chacun des V sous-ensembles est utilisé à tour de rôle comme ensemble de
validation, le reste étant donc utilisé pour l’apprentissage : on calcule, pour
−I
chaque j ∈ {1, . . . , V }, l’ensemble des estimateurs (r̂h j )h∈H fabriqués avec
(Xi , Yi )i∈I
/ j . Ensuite le risque de prédiction pour une fenêtre h est estimé par
V
1 X 1 X
R̂(h) = (Yi − r̂−Ij (Xi ))2
V j=1 nj i∈Ij
où on a noté nj = Card(Ij ).
Dans la pratique on choisit souvent V = 5 ou V = 10.
— Leave-one out : cas particulier du V -fold CV avec V = n.
— Leave-q-out : tout sous-ensemble de cardinal q de l’échantillon est utilisé
comme ensemble de validation et le reste comme ensemble d’apprentissage.
On choisit
ĥ = arg min R̂(h)
h∈H
119
Explicitons un peu plus le cas particulier du "leave-one out". Pour chaque valeur
h de la grille de valeurs H et pour chaque i ∈ {1, . . . , n}, on construit un estimateur
(−i)
r̂h en utilisant toutes les observations sauf la ième. La ième observation est ensuite
(−i) (−i)
utilisée pour mesurer la performance de r̂h par (Yi − r̂h (Xi ))2 . On pose donc
n
1X (−i)
R̂(h) = (Yi − r̂h (Xi ))2 .
n i=1
On minimise R pour trouver ĥ.
Dans la suite on explicite les calculs pour voir le problème de dépendance lié au
risque empirique.
On note X1n = (X1 , . . ., Xn ) et Y1n = (Y1 , . . . , Yn ).
On cherche h tel que E (Y − r̂h (X))2 soit minimal. Remarquez que l’on pourrait
comparer aussi des estimateurs de nature différente. On fait donc disparaitre la
dépendance à h dans la notation.
Si g = r̂, g n’est plus fixe, mais dépend des données (X n , Y n ) et on a
n
1X
E (Yi − r̂(Xi ))2 6= E (Y − r̂(X))2
n i=1
En effet on a, si l’estimateur est symétrique en ses variables (ce qui semble
raisonnable et est le cas des estimateurs par polynômes locaux)
n
1X
E (Yi − r̂(Xi ))2 = E (Y1 − r̂(X1 ))2
n i=1
On indique la dépendance de r̂ à (X n , Y n ) en écrivant r̂(x) = g(X1 , . . . , Xn , Y1 , . . . , Yn , x).
On rappelle qu’on a noté f la densité du couple (X, Y ). On a alors
2 2
E (Y1 − r̂(X1 )) = E (Y1 − g(X1 , . . . , Xn , Y1 , . . . , Yn , X1 ))
Z
= (y1 − g(x1 , . . . , xn , y1 , . . . , yn , x1 ))2 f (x1 , y1 ) . . . f (xn , yn )dx1 dy1 . . . dxn dyn
Tandis que
E (Y − r̂(X))2 = E (Y − g(X1 , . . . , Xn , Y1 , . . . , Yn , X))2
Z
= (y − g(x1 , . . . , xn , y1 , . . . , yn , x))2 f (x1 , y1 ) . . . f (xn , yn )f (x, y)dx1 dy1 . . . dxn dyn dxdy
2
Le risque empirique est un mauvais estimateur du "vrai" risque E (Y − r̂(X)) .
Si (Xn+1 , Yn+1 ) est une nouvelle donnée indépendante de (X1n , Y1n ) et de même
loi que (X, Y ), on a
h i h i
E (Yn+1 − r̂(Xn+1 ))2 = E (Yn+1 − g(X1 , . . . , Xn , Y1 , . . . , Yn , Xn+1 ))2 =
Z
(yn+1 − g(x1 , . . . , xn , y1 , . . . , yn , xn+1 ))2 f (x1 , y1 ) . . . f (xn , yn )f (xn+1 , yn+1 )dx1 dy1 . . . dxn dyn dxn+1 dyn+1
Z
= (y − g(x1 , . . . , xn , y1 , . . . , yn , x))2 f (x1 , y1 ) . . . f (xn , yn )f (x, y)dx1 dy1 . . . dxn dyn dxdy
h i
= E (Y − r̂(X))2
120
On a finalement juste utilisé le fait que
p 2
1X
Yn+k − r̂(Xn+k )
p k=1
particulier sans biais) de E (Y − r̂n−1 (X))2 qui est le "vrai" risque de l’estimateur
r̂n−1 fabriqué à partir de n − 1 données (on s’attend à ce que E (Y − r̂n−1 (X))2
2
soir proche de E (Y − r̂n (X)) où r̂n est l’estimateur de départ, fabriqué avec n
données).
(−i)
On admet la proposition suivante, qui relie les poids associés à l’estimateur r̂h
à ceux associés à l’estimateur r̂h .
Pn (−i) P
Proposition 5.15. Si r̂h (x) = i=1 ωi,h (x)Yi et, pour 1 ≤ i ≤ n, r̂h = j6=i ω̃j,h (x)Yj
alors, pour tout j 6= i
ωj,h (Xi )
ω̃j,h (Xi ) =
1 − ωi,h (Xi )
Remarque 5.16. Cette proposition est également vérifiée pour d’autres estimateurs
linéaires (par exemples les splines).
(−i)
Pour calculer (r̂h )1≤i≤n dans le cas des polynômes locaux, on n’a donc pas
besoin de faire de calculs supplémentaires. Grâce à la proposition précédente on a
facilement le résultat suivant.
121
Pn
Proposition 5.17. Si r̂h (x) = i=1 ωi,h (x)Yi alors
n
1X Yi − r̂h (Xi ) 2
CV (h) =
n i=1 1 − ωi,h (Xi )
Démonstration. On a
n
1X (−i)
CV (h) = (Yi − r̂h (Xi ))2
n i=1
avec
(−i) X
Yi − r̂h (Xi ) = Yi − ω̃j,h (Xi )Yj
j6=i
X ωj,h (Xi )
= Yi − Yj
j6=i 1 − ωi,h (Xi )
P
(1 − ωi,h (Xi ))Yi − j6=i ωj,h (Xi )Yj
1 − ωi,h (Xi )
Pn
Yi − j=1 ωj,h (Xi )Yj
=
1 − ωi,h (Xi )
Yi − r̂h (Xi )
=
1 − ωi,h (Xi )
Il existe une alternative qui consiste à remplacer les ωi,h (xi ) par leur moyenne.
P
Cette alternative s’appelle la validation croisée généralisée. : on pose Ω = ni=1 ωi,h (xi )
puis
n
1X Yi − r̂h (xi ) 2 1 1X n 2
GCV (h) = = Y i − r̂ (x
h i )
n i=1 1 − Ω/n (1 − Ωn )2 n i=1
On minimise ensuite GCV par rapport à h.
Remarquons que si Ω n alors (1 − Ωn )−2 ≈ 1 + 2 Ωn et donc
n 2
1X 2Ω
GCV (h) ≈ Yi − r̂h (xi ) 1 +
n i=1 n
Code R et exemples
On illustre la méthode des polynômes locaux avec une simulation. La fonction
utilisée s’appelle locpoly et appartient au package Kernsmooth. On peut aussi
obtenir une estimation de la fenêtre idéale par la fonction dpill. On va représenter
les résultats associés à diverses fenêtres (une fenêtre sur-lissant, une sous-lissant, et
la fenêtre calculée par la fonction dpill associée à un noyau gaussien). Un noyau
gaussien est utilisé et cette fonction ne permet que l’estimation d’une fonction à une
seule variable. Possibilité d’estimer une dérivée avec l’argument drv (mis à zéro par
défaut) ou bien une densité. Le degré du polynôme correspond à l’argument degree
(par défaut à 1).
Simulation d’un échantillon associé à une fonction r :
122
>x <- seq(0,1,0.05)
>r <- function(x){0.5 + 0.4*sin(2*pi*x)}
>set.seed(10)
>y <- r(x) + rnorm(n=length(x), sd=0.05)
>par(mfrow=c(2,2))
>plot(x, y, pch=16,main="échantillon+ fonction r")
>xtemp <- seq(0,1,0.01)
>lines(xtemp, r(xtemp), lty=2, lwd=2)
Prediction avec la fonction locpoly : on ne peut pas définir une grille de prédiction
quelconque avec cette fonction, seulement une grille de points espacés uniformément
library(KernSmooth)
>h=dpill(x,y) # calcul d’une fenêtre "idéale"
>fenetres=c(0.02,0.25,h)
>for (i in fenetres) {
plot(locpoly(x, y, bandwidth=i,gridsize=101),ylab=paste("h=",i),xlab="",
lwd=2,main="locpoly")
}
123
échantillon+ fonction f locpoly
● ●
●●●
●●
● ●● ●
● ●●
●
0.8
●
0.8
● ● ● ●●●
●●
● ● ●
●
● ●
● ●
● ●
● ● ●
● ●
●● ●●
●
0.6
●
0.6
● ● ● ●●
●
● ●
h= 0.02
● ●
●
● ● ●
● ●
● ● ●
y
● ● ●
● ●
0.4
0.4
● ● ● ●
● ● ●
●
●
● ●
● ● ●
●
● ●
● ●
●● ●
●
●
0.2
0.2
●● ●
●● ●
● ● ●
● ●
●● ●●
● ●● ●
●
●● ●
● ●
●●
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
locpoly locpoly
●●
●●●
●●●
●●●
●● ●
●●●
●●●
●
0.7
●
● ●
●● ●● ●
0.8
●
●● ●
●● ● ●
●
●
●● ●
● ●● ●●
●
h= 0.0719553392537308
●
● ● ●●
●● ●● ●
●● ● ●
●
● ● ●
0.6
● ● ●
●
●
● ● ●
●● ● ●
● ● ●
0.6
●● ● ●
● ● ●
● ● ●
0.5
h= 0.25
●
● ● ●
● ● ●
● ● ●
●● ● ● ●
● ● ●
●● ● ●
● ● ●
0.4
● ● ●
0.4
●● ● ●
●● ● ●
● ● ●
●● ● ●
● ● ●
● ●
0.3
●● ● ●
●● ●● ●
●
● ● ●
●● ● ●
0.2
●● ●
● ●●
●
●● ●● ●
● ●
0.2
●●● ● ●●● ●
●
●●●
●●●
●●●
●●●
●●● ●●●
●●●
●●
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Il existe aussi une fonction appelée loess du package stats qui permet aussi
l’estimation par polynômes locaux, et ce jusqu’à la dimension 4 (de toute façon ce
n’est pas très raisonnable d’aller plus loin en dimension).
Et enfin, il existe aussi le package locfit dont voici quelques paramètres : deg
pour le degré du polynôme local (à 2 par défaut, on utilise rarement au-delà de 3)
, kern pour le noyau (tricube par défaut) , deriv pour estimer une dérivée de la
fonction de régression.
Le choix de la fenêtre est régi par le paramètre alpha. Si on met alpha=c(0,h)
ça donne un estimateur avec une fenêtre égale à h.
Par exemple, si on veut le polynôme local de degré 1 associé à la régression
d’une variable y sur deux variables explicatives x et z avec une fenêtre égale à 0.5,
on utilise locfit(y∼x+z,deg=1,alpha=c(0,0.5)).
Si resultat=locfit(..) alors fitted(resultat) donne les r̂(Xi ) et
residuals(resultat) donne les résidus r(Xi ) − r̂(Xi ).
On va illustrer l’utilisation de la fonction gcvplot associée au package locfit,
fonction qui calcule la validation croisée généralisée pour une série de valeurs de
alpha et fait le graphique correspondant (attention, en abscisse, ce ne sont pas les
valeurs de alpha).
Pour cela on va utiliser les mêmes données simulées.
On utilise une grille de 30 valeurs pour la fenêtre :
>alphamat= matrix(0,ncol=2,nrow=30)
>alphamat[,2]= seq(from=0.1,to=0.8,length=30)
>gcvs= gcvplot(y∼x ,alpha=alphamat,maxk=1000)
La fonction gcvplot est telle que gcvplot$values contient les valeurs de la
validation croisée généralisée (GCV en anglais) et gcvplot$alpha contient les va-
leurs de alpha correspondantes. Donc gcvs$values == min(gcvs$values) donne
la ligne i correspondant à la valeur minimale de la GCV, et avec gcvs$alpha[i,2]
on obtient la valeur de la fenêtre correspondante. Il se peut que plusieurs valeurs
donnent le minimum, auquel cas on prend souvent la plus grande fenêtre donnant
le miminum :
>optband= max(gcvs$alpha[gcvs$values == min(gcvs$values),2])
On peut ensuite fabriquer l’estimateur correspondant à cette fenêtre :
>locfitopt= locfit(y∼x,alpha=c(0,optband),maxk=1000)
>plot(locfitopt,main="locfit fenêtre GCV opt+fonction")
>lines(xtemp,r(xtemp),col=’red’)
125
locfit fenêtre GCVopt+fonction
0.8
0.6
y
0.4
0.2
x
Il y a aussi la possibilité de spécifier une fenêtre différemment, qui n’est pas une
fenêtre constante : pour chaque x où la fonction est évaluée, on utilise une fenêtre
hx telle que qu’il y ait une fraction donnée des Xi dans [x − hx , x + hx ] (ou dans la
boule de centre x et de rayon hx si on est en dimension > 1). Par exemple, si on met
alpha=0.5, on utilise toujours la moitié des données dans l’intervalle [x−hx , x+hx ]).
Ce type de choix est censé être adapté au cas où le design n’est pas distribué assez
uniformément et où on peut avoir peu de données à certains endroits.
127
et donc l’estimateur suivant pour la fonction de régression
n XN
1X
r̂n,N = Yi φj (i/n)φj .
n i=1 j=1
Le choix de la base s’apparente plus au choix du noyau. Les bases les plus fré-
quemment utilisées sont la base trigonométrique et les bases d’ondelettes.
Base Trigonométrique (de Fourier). Elle est donnée par
√ √
φ1 ≡ 1, φ2k : x → 2 cos(2πkx), φ2k+1 : x → 2 sin(2πkx), ∀k ≥ 1.
128
Chapitre 6
Bibliographie conseillée
129
Index
erreur
de première espèce, 13
de seconde espèce, 13
de test, 13
de type I, 13
de type II, 13
hypothèse
composite, 14
simple, 14
test
erreur de, 13
130
Bibliographie
[CHJ+ 12] Pierre-André Cornillon, François Husson, Nicolas Jégou, Eric Matzner-
Lober, and Collectif. Statistiques avec R. PU Rennes, Rennes, 3e édition
revue et augmentée edition, May 2012.
[Dal08] Peter Dalgaard. Introductory Statistics with R. Springer Science & Busi-
ness Media, August 2008.
[Loa99] Clive Loader. Local Regression and Likelihood. Springer, New York, 1999
edition edition, July 1999.
131