MATH 444 – Statistique multivariée

Dr. Jean-Marie Helbling

Notes de cours par

Léo Belzile
[Link]@[Link]

V ERSION DU 2 MAI 2018

H IVER 2015 ET 2016, EPFL
Écrire un courriel à l’auteur si vous trouvez une coquille.
Ces notes n’ont été que partiellement révisées et devraient être consultées avec prudence.
Licencié sous Creative Commons Attribution-Non Commercial-ShareAlike 3.0 Unported
 Table des matières
1 Introduction 4

1.1 Jeux de données . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2 Organisation des données . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Rappel de probabilités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Loi normale multivariée et ses dérivées 8

2.1 La loi normale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Estimation des paramètres µ et Σ . . . . . . . . . . . . . . . . . . . . . . . 12

2.3 Lois dérivées de la loi normale multivariée . . . . . . . . . . . . . . . . . 13

2.4 Vérification de la normalité multivariée . . . . . . . . . . . . . . . . . . . 17

3 Inférences relatives aux paramètres de lois normales 23

3.1 Distribution de X̄ et S . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.2 Tests relatifs à une moyenne . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.3 Tests relatifs à une variance . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.4 Comparaison de deux ou plusieurs populations normales . . . . . . . . . 30

4 Analyse en composantes principales 36

4.1 Inférence en ACP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4.2 Interprétation et qualité d’une ACP . . . . . . . . . . . . . . . . . . . . . . 40

4.3 ACP, décomposition en valeurs singulières et régression . . . . . . . . . 42

5 Analyse canonique 49

5.1 Définition et propriétés de l’analyse canonique . . . . . . . . . . . . . . . 49

5.2 Inférence en analyse canonique . . . . . . . . . . . . . . . . . . . . . . . . 51

5.3 Généralisation : analyse canonique et tableaux de contingence . . . . . . 52

6 Partitionnement de données 57

6.1 Méthodes de regroupement hiérarchiques . . . . . . . . . . . . . . . . . . 57

6.2 Méthodes non hiérarchiques . . . . . . . . . . . . . . . . . . . . . . . . . . 59

2
 6.3 Méthodes basées sur un modèle . . . . . . . . . . . . . . . . . . . . . . . . 63

6.4 Algorithme d’espérance-maximisation . . . . . . . . . . . . . . . . . . . . 64

7 Analyse discriminante et classification 73

7.1 Le problème général de la classification . . . . . . . . . . . . . . . . . . . 73

7.2 Analyse discriminante et classification par la méthode de Fisher . . . . . 76

7.3 Autres considérations sur la discrimination . . . . . . . . . . . . . . . . . 80

7.4 Évaluation de la qualité de la discrimination . . . . . . . . . . . . . . . . 80

7.5 Discrimination et régression logistique . . . . . . . . . . . . . . . . . . . . 81

7.6 Arbres de décisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

8 Copules 86

8.1 Motivation et définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

8.2 Théorème de Sklar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

8.3 Mesure de la dépendance . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

8.4 Familles et modèles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

8.5 Estimation et inférence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

3
 Chapitre 1
Introduction
Pour motiver le cours, on regarde quelques jeux de données qui seront utilisés dans le
cours.

1.1 Jeux de données

• Crabes :
– distribution de données
– test de différences de moyennes (même longueur de carapace chez les mâles et les
femelles).
– séparer au mieux les groupes pour classer de nouveaux cas (MANOVA).
• Examens
– analyse canonique : liens entre les types d’examen (livre ouvert ou fermé, par
exemple)
– représentation la plus simple possible ; par régression, ou projection
– représentation en dimension inférieure (analyse en composante principale)
• Composition chimique du sang (BMDP)
– discrimination
– régression (dépendance du cholestérol et de l’âge)
– données manquantes
– aberrances
• Ventes
– lien entre les notes d’examen et les ventes
– formation de groupes plus ou moins homogènes (agglomération, ou “clustering”)

1.2 Organisation des données

Les données seront stockées sous forme matricielle de dimension n × p. On considère

un échantillon de taille n sur lequel on mesure p variables.

X11 X1p
 
···
 .. 
X . X2p 
X =  .21
 
 . .. .. 
 . . . 


Xn1 ··· Xnp

et les observations sont indépendantes et identiquement distribuées (iid). On dénote

cette matrice X = ( X 1> , X 2> , . . . , X > > e
n ) où le vecteur colonne X i est la i observation.
Sur ce tableau de données, on calcule les statistiques descriptives suivantes :

• la moyenne X̄ = ( X̄1 , X̄2 , . . . , X̄ p )> = n−1 X> 1n , où X̄1 correspond à la moyenne de

la 1e variable sur tous les individus (moyenne de la 1e colonne de X).

4
 • les variances :
n
1
S2j = ∑ ( X − X̄ j )2 , j = 1, . . . , p
n − 1 i=1 ij

la variance de la variable j sur tous les individus.

• les covariances
n
1
n − 1 i∑
S jk = ( Xij − X̄ j )( Xik − X̄k ), j 6= k
=1

ce qui donne la matrice de variance-covariance

S12 S jk S1p
 
··· ···
 . .. .. .. 
 .. . . ··· . 
S= .
 
 . .. .. 
 . Skj . ··· . 


S p1 ··· ··· ··· S2p

n
1
n − 1 i∑
= ( X i − X̄ )( X i − X̄ )> .
=1
1  > 
= X X − n X̄ X̄ >
n−1
1 1 > >
 
>
= X X − X 11 X
n−1 n
1 1
 
> >
= X I − 11 X
n−1 n

Noter que I − 11> /n est une matrice idempotente, et de ce fait son rang est égal à sa
trace. On a donc

1 1
 
tr I − 11> = n − n = n − 1.
n n
q
• la matrice des corrélations, avec élément hors-diagonal r jk = S jk / S2j Sk2 .

 
1 ··· r1p
 . .. .. 
 ..
R= . . 

r p1 ··· 1
1 1
= (diag (S))− 2 S (diag (S))− 2

• notions de variance généralisée : det(S) ou det( R), tr(S) ou valeurs propres.

Note
Les valeurs propres de S sont positives puisque S est symétrique (et définie positive).

5
1.3 Rappel de probabilités

Soit X = ( X1 , . . . , X p )> , un vecteur aléatoire de dimension p, avec µ = E ( X ) =

(E ( X1 ) , . . . , E X p )> et

   
Σ X = Cov ( X ) = E ( X − µ)( X − µ)> = E X X > − µµ> .

Proposition 1.1 (Transformations linéaires)

Si Z = CX et W = DX, alors

E ( Z ) = µ Z = Cµ X ,
Cov ( Z ) = Σ Z = CΣ X C>
et
Cov ( Z, W ) = CΣ X D> .

Par exemple, Cov a> X, b> X = a> Σ X b.

Définition 1.2 (Fonction de répartition (multivariée))

La fonction de répartition F : Rd → [0, 1] est définie pour tout x1 , . . . , xd ∈ R par

F ( X1 , . . . , X n ) = P ( X1 ≤ x 1 , X2 ≤ x 2 , . . . , X d ≤ x d ) .

Propriétés de F
F caractérise le comportement aléatoire de x. Toute fonction de répartition F satisfait
les propriétés suivantes.
(i)

lim . . . lim F ( x1 , . . . , xd ) = 1;
x1 → ∞ xd →∞

(ii)

lim F ( x1 , . . . , xd ) = 0, ∀ xi
xi →−∞

puisque l’intersection de ∅ avec n’importe lequel autre ensemble est ∅, de probabilité

nulle.
(iii) F est non-décroissante et continue à droite dans chacun de ces d arguments. Par
example, on peut considérer fixer x2 , . . . , xd ∈ R, de telle sorte que t 7→ F (t, x2 , . . . , xd )
est non-décroissante et continue à droite.

6
(iv) F est d−monotone. Dans le cas d = 2, cela signifie que pour x1 < x1∗ et x2 < x2∗

F ( x1∗ , x2∗ ) − F ( x1 , x2∗ ) − F ( x1∗ , x2 ) + F ( x1 , x2 ) ≥ 0

F IGURE 1 – 2−monotonicité illustré par le biais d’ensemble

b2

B D

a2

A C

a1 b1

L’image illustre P( A ∪ B ∪ C ∪ D ) − P( A ∪ C ) − P( A ∪ B) + P( A) puisqu’il s’agit d’ensembles

disjoints. En terme d’aire, c’est ( A + B + C + D ) − ( A + C ) − ( A + B) + A = D et ainsi la proba-
bilité P( a1 < X1 < b1 , a2 < X2 < b2 ) = D est positive.

Le cas d = 3 donne

F ( x1∗ , x2∗ , x3∗ ) − F ( x1 , x2∗ , x3∗ )+ F ( x1 , x2 , x3∗ ) − F ( x1 , x2 , x3 )

− F ( x1∗ , x2 , x3∗ )+ F ( x1 , x2∗ , x3 )
− F ( x1∗ , x2∗ , x3 )+ F ( x1∗ , x2 , x3 ) ≥ 0

et en dimension d arbitraire

∑ (−1)v(c) × F (c1 , . . . , cd ) ≥ 0
ci ∈{ xi ,xi∗ }
i ∈{1,...,d}

où v(c) est le nombre de constantes ci sans étoiles. Inversement, toute fonction F :

Rd → [0, 1] qui satisfait (i)–(iv) est une fonction de répartition valide, c’est-à-dire qu’il
existe une mesure de probabilité qui engendre F. La preuve est la même que dans le
cas univarié ; l’élément clé est que ( x1 , x1∗ ] × ( x2 , x2∗ ] × · · · × ( xd , xd∗ ] sont préciséments
des générateurs de la tribu borélienne sur Rd .

Les propriétés (i)–(iii) sont suffisantes en une dimension, mais la propriété essentielle
pour la construction d’une fonction de répartition multivariée est la (iv). On dénomme
loi marginale la loi correspondant à une sous-composante de dimension d − q du d-
vecteur X, pour 1 ≤ q ≤ d.

7
 Chapitre 2
Loi normale multivariée et ses dérivées
2.1 La loi normale
Définition 2.1 (Loi normale univariée)
On caractérise la loi normale par sa densité

2 !
1 1 x−µ

f (x) = √ exp −
2πσ2 2 σ

où ( x − µ)2 /σ2 = ( x − µ)(σ2 )−1 ( x − µ). On dénote par X ∼ N (µ, σ2 ).

Définition 2.2 (Loi normale multivariée)

La densité de la loi multinormale est

1 1
 
f ( x) = 1
exp − ( x − µ)> Σ−1 ( x − µ)
(2π ) p/2 |Σ| 2 2

Note
La loi multinormale est entièrement caractérisée par ses deux premiers moments, res-
pectivement µ et Σ. Son utilisation pour certaines quantités d’intérêt (la moyenne arith-
métique) est justifiée asymptotiquement par le théorème central limite.

Proposition 2.3
Si Σ est positive définie, alors E ( X ) = µ et Var ( X ) = Σ.

Preuve Puisque Σ est positive définie, il existe une matrice orthogonale P telle que
P> ΣP = Λ, où Λ = diag (λ1 , . . . , λ p ) est une matrice diagonale de valeurs propres.
Posons X = PY + µ et calculons la densité de Y. Le jacobien de la transformation est
p
J = ∂X i /∂y j = |P| et l’on note que |Σ| = |P|2 ∏i=1 λi . La densité de Y est égale à
i,j
f Y (y) = f X (Py + µ)|J|. Dès lors,

1 1
 
f (y) = p/2
exp − y> P> Σ−1 Py
(2π ) 2
p
λ1 · · · λ p
1 1 > −1
 
= exp − y Λ y
(2π ) p/2 λ1 · · · λ p 2
p

= f ( y1 ) f ( y2 ) · · · f ( y p )

puisque Σ−1 = (PΛP> )−1 = PΛ−1 P> . Ainsi, Y ∼ N p (0, Λ) et en plus, Yi ∼ N (0, λi ).
Par linéarité, E ( X ) = PE (Y ) + µ = µ. L’on a également

Var ( X ) = Var (PY + µ) = Var (PY ) = PVar (Y ) P> = PΛP> = Σ.

8
 

Pour la loi normale multivariée, les courbes de densité constante sont définies par c2 =
√
( X − µ)> Σ−1 ( X − µ). Elles représentent des ellipses de centre µ et d’axes ±c λi vi où
λi est une valeur propre de Σ avec vecteur propre vi . La forme dépend de Σ et c. 1

F IGURE 2 – Ellipses de confiance pour la distribution marginale {analyse, statistique} du

jeu de données examens

Ellipses de confiances
100

0.95
80

0.85

0.65
0.45
60
Analyse

0.25

0.05
40

0.15
0.35
0.55
0.75
20
0

0 20 40 60 80 100

Statistique

Proposition 2.4
La forme quadratique

Q = ( X − µ)> Σ−1 ( X − µ) ∼ χ2p

suit une loi de khi-carré avec p degrés de liberté.

Définition 2.5 (Fonction caractéristique)
On définit la fonction caractéristique d’un vecteur aléatoire comme
  Z
ΦX (t) = E exp(it X ) = exp(it> x) f X ( x) dx
: >

1. Il en résulte qu’il est facile de calculer la probabilité pour une loi multinormale d’être dans une ellip-
soïde, contrairement par exemple à un parallélépipède.

9
pour tout t ∈ R p . Ses propriétés sont

• Existence : la fonction caractéristique existe toujours et elle satisfait ΦX (0) = 1 et

|ΦX (t)| ≤ 1.
• Unicité : deux vecteurs aléatoires dotés de la même fonction caractéristiques si et
seulement si ils ont la même distribution.
• Inversion : si ΦX (t) est absolument intégrable, alors la densité de X est donnée par

1
Z
f (x) = exp(−it> x)ΦX (t) dt
(2π ) p

• Indépendance : la fonction caractéristique du vecteur X = ( X 1 X 2 ) est telle que

ΦX (t) = ΦX 1 (t)ΦX 2 (t) si et seulement si X 1 ⊥⊥ X 2.
• Moments :  
j1 +···+ j p

j jp
 ∂
E X11 X p = i−( j1 +···+ j p )  ΦX ( t ) 
j1 j p
∂t1 ∂t p
t =0

• La fonction caractéristique du (d − q)-vecteur marginal issu d’un d-vecteur X pour

t
1 ≤ q ≤ d est ΦX ( 01q )
• Si X et Y sont indépendants, alors ΦX +Y (t) = ΦX (t)ΦY (t)

Théorème 2.6 (Cramér-Wold)

La distribution d’un vecteur aléatoire X est complètement déterminée par l’ensemble
de toutes les lois univariées formées des combinaisons linéaires a> X, pour tout a ∈ R p .

Corollaire 2.7
X suit une loi multinormale si a> X est (multi)normale pour tout a ∈ R p .

Proposition 2.8 (Fonction caractéristique de la loi multinormale)

La fonction caractéristique de la loi multinormale est donnée par

1
 
exp it> µ − t> Σt
2

Preuve On procède au calcul en complétant le carré :

 
ΦX (t) = E exp(it> X )
1
Z    
− p/2 −1/2 > −1
= (2π ) |Σ| exp − ( x − µ) Σ ( x − µ) exp it> x dx
2
1 > −1
 
y= x−µ
Z  
− p/2 −1/2
= (2π ) |Σ| exp − y Σ y exp it> (y + µ) dy
2
1
Z  
− p/2 −1/2 > −1
= (2π ) |Σ| exp − (y − iΣt) Σ (y − iΣt) dy
2
1
 
× exp it> µ + t> Σti2
2

10
puisque

1 1 
− y> Σ−1 y + it> y + it> µ = − y> Σ−1 y − 2it> ΣΣ−1 y + it> µ
2 2
1 1
= − (y − iΣt)> Σ−1 (y − iΣt) + it> µ + t> Σti.
2 2


Corollaire 2.9
Si X ∼ N p (µ, Σ) et qu’on partitionne le vecteur en X 1 , X 2 , alors
! ! !!
X1 µ1 Σ11 Σ12
X= ∼ Np ,
X2 µ2 Σ21 Σ22

Note
En supposant que X suit une loi normale, les sous-vecteurs X 1 , X 2 de dimension res-
pective q, d − q sont indépendants 2 (dénoté X 1 ⊥ > =
⊥ X 2 ) si et seulement si Σ12 = Σ21
Oq,d−q .
Proposition 2.10 (Distribution de la moyenne)
La distribution de la moyenne arithmétique de multinormales est X̄ ∼ N p (µ, n−1 Σ).

iid
Preuve On écrit X̄ = n−1 ∑in=1 X i avec X i ∼ N p (µ, Σ). Puisque X̄ est une combinaison
linéaire des X i , on a par linéarité E ( X̄ ) = µ et Var ( X̄ ) = n−2 nΣ = n−1 Σ. 

Proposition 2.11
Si X ∼ N p (µ, Σ) et Y = CX + d, où C est une matrice non-stochastique q × p et d un
vecteur colonne de dimension q. Alors Y ∼ Nq Cµ + d, CΣC>

Corollaire
 2.12
I O

Si C = Oq−s s, s Oq−s,s,q−q−s s pour 1 ≤ s ≤ q et d = 0q , alors on obtient les distributions
marginales, lesquelles sont normales. Si on partitionne X = ( X 1 , X 2 )> de dimensions
 I Os,p−s
s

s, p − s respectivement, alors l’utilisation de la matrice C = −Σ Σ−1 I dans le
21 11 p−s
précédent résultat donne
! ! !!
X1 µ1 Σ11 Os,p−s
Y≡ −1 ∼ Np −1 , −1
X 2 − Σ21 Σ11 X1 µ2 − Σ21 Σ11 µ1 O p−s,s Σ22 − Σ21 Σ11 Σ12

−1
et donc la distribution conditionnelle est X 2 | X 1 ∼ N p−s µ2 + Σ21 Σ11 ( X 1 − µ1 ), Σ22.1 ,


−1 −1
où Σ22.1 := Σ22 − Σ21 Σ11 Σ12 . En effet, sachant X 1 , −Σ21 Σ11 X 1 est une constante et on
−1
peut déduire E ( X 2 ) et Var ( X 2 ) de l’expression pour X 2 − Σ21 Σ11 X 1.
2. La normalité des marges n’est pas une condition suffisante ; la preuve se fait par le biais de la fonction
caractéristique.

11
2.2 Estimation des paramètres µ et Σ

Dans cette section, on dérive le maximum de vraisemblance des paramètres de la loi

multinormale. Soit un échantillon X 1 , . . . , X n . La log-vraisemblance ` = log( L) est don-
née par

np n   1 n
`( X; µ, Σ) = − log(2π ) + log det Σ−1 − ∑ ( X i − µ)> Σ−1 ( X i − µ)
2 2 2 i =1
n
!
np n   1
=− log(2π ) + log det Σ−1 − tr Σ−1 ∑ ( X i − µ)( X i − µ)>
2 2 2 i =1

On s’intéresse à la trace et tout particulièrement à la somme des termes centrés. L’on a

pour cette dernière

n n
∑ (X i ± X̄ − µ)(X i ± X̄ − µ)> = ∑ (X i − X̄ )(X i − X̄ )> + n(X̄ − µ)(X̄ − µ)>
i =1 i =1

puisque les termes croisés ∑in=1 ( X i − X̄ )( X̄ − µ)> sont nuls. La log-vraisemblance par
rapport à µ est proportionnelle à

n
!!
µ 1
`( X; µ, Σ) ∝ − tr Σ−1 ∑ (X i − X̄ )(X i − X̄ ) >
+ n( X̄ − µ)( X̄ − µ) >
2 i =1

et en simplifiant l’écriture, on a que

np n  
`( X; µ, Σ) = − log(2π ) + log det Σ−1
2 2
1  −1  n
− tr Σ S XX − ( X̄ − µ)> Σ−1 ( X̄ − µ).
2 2
Le maximum de vraisemblance de la moyenne découle immédiatement de la dernière
expression et on trouve µ
b = X̄. En dérivant par rapport à la matrice de précision, on
obtient
∂` n 1
= Σ − S>
∂Σ−1 2 2 XX

b = n−1 S XX . On procède maintenant au calcul du biais

et en égalant à zéro, on obtient Σ
des estimateurs du maximum de vraisemblance (EMV) :

n
!
1 n
b) = E
E (µ
n ∑ Xi =
n
µ=µ
i =1
n
!
b = 1E
 
E Σ
n ∑ (X i − X̄ )(X i − X̄ ) >
i =1

12
 n
!
1
= E
n ∑ (X i − µ)(X i − µ) >
− n(µ − X̄ )(µ − X̄ ) >
i =1
1
= nΣ − Var ( X̄ )
n
n−1
= Σ
n

et ce dernier estimateur a un biais de n−1 Σ. On peut prendre

n
1
S := ∑ ( X − X̄ )( X i − X̄ )>
n − 1 i =1 i

comme estimateur de Σ.

2.3 Lois dérivées de la loi normale multivariée

Définition 2.13 (Loi de Wishart)
La distribution de Wishart, dénotée W est une distribution sur les matrices aléatoires
positives définies et une généralisation multivariée de la distribution de χ2 .

On considère n vecteurs aléatoires indépendants X 1 , . . . , X n avec X i ∼ N p (µi , Σ). Alors

n
W= ∑ X i X i> = X> X ∼ W p (n, Σ, M )
i =1

où M = (µ1> , . . . , µ> >

n ) est une matrice p × n. La loi est dite centrée si M = O et on
écrit alors W p (n, Σ).

Proposition 2.14 (Propriétés de la loi de Wishart)

1. E (W ) = nΣ + M M >
·
2. Si c est un vecteur non aléatoire, alors c> W c ∼ σc2 χ2n avec σc2 = c> Σc. Si C est une
matrice (non stochastique) de rang q et que W ∼ W p (n, Σ), alors

CW C> ∼ Wq (n, CΣC> ).

3. X > AX ∼ W p (rang(A), Σ, ·) si et seulement si A2 = A. De plus, la distribution W p

est centrée si AM = 0.
4. X > AX et X > BX de lois Wishart sont indépendants si et seulement si AB = 0
5. X > B et X > AX de lois Wishart sont indépendants si et seulement si B> A = 0.
6. Si W ∼ W p (n, Σ), alors |W |/|Σ| ∼ χ2n χ2n−1 · · · χ2n− p+1 où les p variables aléatoires
khi-carrées sont indépendantes.
7. Si W ∼ W p (n, Σ) alors pour tout vecteur a fixé, on a

a > Σ −1 a
∼ χ2n− p+1 .
a > W −1 a

13
 on obtient 1> W1 ∼ χ2n− p+1 .

Preuve

1. Écrire
n n
W= ∑ X i X i> = ∑ (X i − µi + µi )(X i − µi + µi )>
i =1 i =1
n n n n
= ∑ (X i − µi )(X i − µi )> + ∑ µi (X i − µi )> + ∑ (X i − µi )µi> + ∑ µi µi>
i =1 i =1 i =1 i =1

et il suffit maintenant de calculer l’espérance. On obtient

n   n
E (W ) = ∑E ( X i − µi )( X i − µi )> + ∑ µi µi>
i =1 i =1
>
= nΣ + M M

2. On a W = ∑in=1 X i X i> avec X i ∼ N p (µi , Σ). En multipliant par une matrice C, on

obtient la forme quadratique

n n
CW C> = ∑ CX i (CX i )> = ∑ Zi Zi> .
i =1 i =1

Or, Z i ∼ Nq Cµi , CΣC> , ce qui implique que CW C> ∼ Wq (n, CΣC> , ·).

Dans le cas vectoriel,

n n n 2
Yi

c> W c = ∑ (c> X i )(X i> c) = ∑ Yi2 = σc2 ∑ σc
= σc2 χ2n
i =1 i =1 i =1

avec Yi ∼ N c> µi , σc2 où σc2 := c> Σc.

3. On prouve un énoncé indirect : posons Z = Xc pour un vecteur constant c 6= 0.

Alors X > AX ∼ W p (r, Σ, ·) avec r = rang(A) = tr(A) si et seulement si Z > AZ ∼
σc2 χ2r , pour tout c.

La nécessité découle de la proposition (point 2). On a X > AX ∼ W p (r, Σ, ·), ce qui

implique que Z > AZ = c> X > AXc ∼ σc2 χ2r . Pour la suffisance, on a Z > AZ/σc2 ∼ χ2r
et Z ∼ N (·, σc2 I p ). On a A est idempotente de rang r que A = ∑ri=1 ai ai> avec ai un
vecteur propre orthonormé de A. Ainsi,

r r
X > AX = ∑ X > ai ai> X = ∑ U i U i> .
i =1 i =1

Puisque les ai sont orthonormaux, YU i ∼ N p (ai> M, Σ) sont indépendants, ce qui

implique que X > AX suit une loi Wishart W p (r, Σ, ·) qui est centrée si AM = 0.

14
 4. Nécessité de la condition : X > AX ∼ W p (r, Σ, ·) et X > BX ∼ W p (s, Σ, ·) sont indé-
pendantes. Alors Z > AZ et Z > BZ sont indépendantes (fonctions de deux lois indé-
pendantes). Par point 3, on a Z > AZ/σc2 ∼ χ2r et Z > BZ/σc2 ∼ χ2s indépendants pour c
quelconque, donc AB = 0.

Suffisance de la condition : par point 3, A et B sont idempotentes de rang r et s, res-

pectivement. On a donc A = ∑ri=1 ai ai> et B = ∑sj=1 b j où ai et b j sont des vec-
teurs propres. De plus, ai> b j = ai> A> b j = 0 pour tout i, j par hypothèse. On a donc
X > AX = ∑ U i> U i avec X > ai et X > BX = ∑ V j V > >
j avec V j = X b j indépendants.

5. Pour la nécessité de la condition, on suppose que X > B ∼ N p (·, ·) et X > AX ∼

W p (n, Σ, ·) sont indépendants. En utilisant le même raisonnement que dans la preuve
précédente, on a Z > B = c> X > c ∼ N 0, σc2 B> B et Z > AZ ∼ σc2 χ2r . De plus, Z > B et

Z > AZ sont indépendants comme fonctions de vecteurs aléatoires indépendants.

Pour la suffisance, on a Z > B ∼ N p (·, ·) et Z > AZ ∼ W p (n, Σ, ·) sont indépendants

pour c quelconque implique X > B ∼ N p (·, ·) et X > AX ∼ W p (n, Σ, ·). Comme Z ∼
N (·, σc2 I ) on a Z > B ∼ N (·, σc2 B> B) et Z > BB> Z/(B> B ∼ σc2 χ21 avec BB> A = 0.
Multiplier par B> , on obtient que b> A = 0 puisque B> B est constante. B est perpen-
diculaire à tous les A implique que U 1 , . . . , U r tels que définis précédemment et X > B
sont indépendants.
6. Par induction. Le cas p = 1 est vrai car W ∼ σ2 χ2n . Supposons que l’énoncé est vrai
pour le cas p ∈ N et que W ∼ W p+1 (n, Σ, ·). On écrit W comme

−1
! ! !
W 11 W 12 W 11 0p Ip W 11 W 12
W= = −1
W 21 W 22 W 21 1 0>
p W 22 − W 21 W 11 W 12

pour W 22 une matrice 1 × 1 et par hypothèse |W 11 | ∼ |Σ11 |χ2n · · · χ2n− p+1 . Or W 22.1 =
−1
W 22 − W 21 W 11 W 12 ∼ W1 (n − p, Σ22.1 ) et indépendant de (W 11 , W 12 ). La décompo-
sition |W | = |W 11 ||W 22.1 | et |Σ| = |Σ11 ||Σ22.1 | assure le résultat.

Définition 2.15 (Loi T 2 de Hotelling)

On fait le parallèle avec la loi normale univariée : si X ∼ N (µ, σ2 ) et W ∼ σ2 χ2n sont
des variables aléatoires indépendantes, alors

X −µ
T = qσ ∼ Tn ,
W
nσ2

soit une loi de Student-t. En général, on a la statistique (n − 1)S2 /σ2 ∼ χ2n−1 pour un

15
échantillon, ce qui implique que

X̄ −µ
 
√
X̄ − µ σ/ n
t= √ = r ∼ T n −1
S/ n S2 ( n −1)
σ 2 ( n −1)

et T 2 = n( X − µ)2 /W a une loi qui est proportionnelle à une loi de Fisher F (1, n).

La loi T 2 de Hotelling est une généralisation multivariée de cette statistique.

T 2 = n( X − µ)> W −1 ( X − µ) ∼ T p,n
2

2 la distribution.
avec X ∼ N p (µ, Σ) et W ∼ W p (n, Σ) indépendants. On dénote par T p,n

Proposition 2.16
La loi de Hotelling T 2 ∼ T p,n
2 est équivalente à une loi de Fisher, soit

n − p + 1 T2
∼ F p,n− p+1
p n

Preuve Esquisse : poser Y = X − µ, ainsi

T2 Y > Σ −1 Y G
= Y > W −1 Y = > −1 ≡
n Y Σ Y H
Y > W −1 Y

où G et H correspondent au numérateur et au dénominateur, respectivement. Posons

maintenant Y = a fixé. La loi conditionnelle de H étant donné Y = a est χ2n− p+1 par
la propriété 7 de la loi de Wishart et donc ne dépend pas de Y (car la distribution ne
dépend pas de a). Ainsi H est indépendant de Y et donc de G (car Σ est connu). La
forme quadratique G = Y > Σ−1 Y suit une loi χ2p et le ratio G/H est le ratio de lois χ2p
sur χ2n− p+1 . On obtient donc (n − p + 1) G/( pH ) ∼ F p,n− p+1 . 

Remarque (Lois elliptiques)

On peut généraliser les résultats obtenus avec la loi normale à une classe de vecteurs
aléatoires dont la densité est de la forme générale f ( x) ∝ g(( x − µ)> Σ−1 ( x − µ)). Les
principaux exemples sont la loi multinormale et la loi de Student-t.
Définition 2.17 (Lois elliptiques)
d
X suit une loi elliptique s’il admet la représentation stochastique X = µ + RAΘ, où

• µ est un vecteur de localisation,

• R est une variable aléatoire positive,
• AA> est la décomposition de Cholesky de la matrice de variance Σ et
• Θ ∼ Sd est une variable uniforme sur la d − 1-sphère.

16
Les variables elliptiques sont dotées de symétrie radiale, qui associe la fonction de ré-
partition centrée à sa fonction de survie. La covariance nulle entre deux marges im-
plique leur indépendance seulement si la loi est multinormale, d’où sa prépondérance
dans les modèles graphiques. Notez que selon la forme du générateur g, il est possible
que certains moments de la distribution n’existent pas. On peut dériver des propriétés
d
à partir de la forme quadratique R2 = ( X − µ)> Σ−1 ( X − µ).

Tableau 1 – Quelques lois elliptiques

Copule FR2 g(t)

1
Gaussienne χ2 ( d ) cd e− 2 u
Student-t F (d, ν)/d cd (1 + t/ν)−(d+ν)/2
IIe fn de Pearson B(d/2, ν + 1) cd (1 − u)ν 1ν>−1,−1≤u≤1
VIIe fn de Pearson m · BII (d/2, b − d/2) cd (1 + t/m)−b 1b>d/2,m>0
Kotz symétrique Γ( a + d/2 − 1, ν) cd ub−1 e−νu 1ν>0,2a−d>2
Laplace cd e−|u|
Logistique c d e − u / (1 + e − u )2

2.4 Vérification de la normalité multivariée

Cette dernière implique que les marges sont normales. On se penche d’abord sur les
tests de normalité univariée.

2.4.1. Normalité univariée

Il existe des méthodes informelles (diagrammes quantiles-quantiles ou diagrammes

probabilité-probabilité) et des tests formels. Au lieu d’effectuer un test formel d’adé-
quation, on peut procéder à une vérification informelle à l’aide de graphiques appelés
diagrammes Q–Q (quantiles–quantiles) et P–P (probabilités–probabilités). Posons F la
fonction de répartition correspondant à la density f
Définition 2.18 (Diagramme quantile–quantile)
Un graphique des rangs des observations, x[i] versus ai avec ai = F −1 ( pi ) et pi = (i −
c)/(n − 2c + 1), un estimé empirique de la probabilité cumulative des observations.
Typiquement, on choisit c = 0 (pseudo-observations) ou c = 1/3 qui donne une fonc-
tion empirique approximativement sans biais médian, puisque alors pi ≈ med( F ( x[i] )).
C’est un paramètre tel que 0 ≤ c ≤ 1, et F −1 représente la fonction inverse de F. Si la
représentation suit une droite, les données sont susceptibles de provenir de F. Noter
que le diagramme Q–Q magnifie les observations extrêmes. L’abscisse peut être stan-
dardisé (dans le cas d’un diagramme Q–Q normal, on parle alors de l’échelle N (0, 1))
ou pas.
Définition 2.19 (Diagramme probabilités–probabilités)
Une représentation graphique de z[i] = F (( x[i] − µ
b)/b
σ ) versus pi tel que défini ci-

17
dessus. Si la représentation est une droite à 45deg , on peut raisonnablement penser que
les données proviennent de F.

Exemple 2.1
On considère un jeu de données du gain de poids (masse) de femelles rats sous un
régime à forte teneur protéique.

Tableau 2 – Données de l’expérience et positions des points des diagrammes Q–Q et

P–P normaux standardisés

xi − µ
 
i xi pi = n+i 1 Φ −1 ( p i )b
σ+µ Φ
b
σ
b b

1 83 0.077 89.50 0.0418

2 97 0.154 98.18 0.1141
3 104 0.231 104.25 0.2272
4 107 0.308 109.25 0.2717
5 113 0.385 113.73 0.3717
6 119 0.462 117.93 0.4813
7 123 0.538 122.07 0.5558
8 124 0.615 126.27 0.5742
9 129 0.692 130.75 0.6630
10 134 0.769 135.75 0.7436
11 146 0.846 141.82 0.8879
12 161 0.923 150.50 0.9724

Diagramme quantile-quantile normalisé Diagramme probabilité-probabilité

1.0
160

0.8
140
Échantillon observé

Φ (( xi − x̄ )/σ̂)

0.6
120

0.4
100

0.2
80

90 100 110 120 130 140 150 0.2 0.4 0.6 0.8

Φ−1 (σ̂ · i/(n + 1) + x̄ ) pi = i/(n + 1)

Proposition 2.20 (Test du χ2 d’adéquation)

Le test d’adéquation du χ2 est basé sur H0 : X ∼ Fθ. Soit X un échantillon n × 1 d’ob-
servations indépendantes et ν le nombre de paramètres. On peut partitionner l’espace
en k intervalles disjoints A1 , . . . , Ak ; dénotons par e j la fréquence théorique de X dans
Ai et o j la fréquence observée. Cela revient à tracer un histogramme dans lequel o j est

18
le nombre d’observations dans la classe j. On a

k ( o j − e j )2 d 2
X2 = ∑ −
→ χ k − ν −1
j =1
ej

quand n → ∞. L’approximation est bonne pour n large (au bas mot 50) et e j > 5. Dans

Histogramme
40
30
Fréquence

20
10
0

-3 -2 -1 0 1 2 3

le cas de la loi normale, on doit calculer µ et σ, ce qui correspond au cas ν = 2.

Si les décomptes e j sont trop petits, il convient de former de plus gros groupes. Si on
dessine l’histogramme, porter attention à faire toujours tous les intervalles de longueur
égale. Le nombre d’observations pour l’application du test est important, et il est pré-
férable de faire un diagramme Q–Q dans les cas où la taille de l’échantillon est faible.

Proposition 2.21 (Test d’adéquation de Kolmogorov-Smirnov)

Comparaison de la fonction de répartition théorique et empirique. On considère la plus
grand distance, soit

sup | Fn ( x ) − F ( x )| .
x

Il est important de noter que tous les paramètres doivent être spécifiés, soit exactement
connus.

Il existe aussi des tests spécifiquement construits pour la loi normale contrairement aux
deux précédents qui sont des tests omnibus, comme les statistiques Shapiro–Wilks et
Anderson–Darling A2n . On donne en exemple la première :

19
 Kolmogorov-Smirnov

1.0
0.8
0.6
Fn ( x )

0.4
0.2
0.0

-2 -1 0 1 2

Proposition 2.22 (Test de Shapiro–Wilks)

Le test de Shapiro–Wilks a comme statistique
2
(n)

∑in=1 an−i+1 X[n−i+1] − X[i]
W=
∑in=1 ( X[i] − X̄ )2

où X[i] est la ie statistique d’ordre (en ordre croissant). Par exemple, X[1] dénote la plus
petite observation. Les valeurs critiques tout comme la loi de W, sont tabulées.

2.4.2. Normalité multivariée

Proposition 2.23 (Distance de Mahalanobis)

On rappelle que ( X − µ)> Σ−1 ( X − µ) ∼ χ2p si X ∼ N p (µ, Σ). On peut utiliser pour n
grand l’approximation

·
d2i = ( X i − X̄ )> S−1 ( X i − X̄ ) ∼ χ2p ,

soit la distance de Mahalanobis entre X i et sa moyenne. Le diagramme Q–Q des d2i

versus les quantiles de la loi χ2 est un test graphique de multinormalité. On n’a en
revanche pas de données indépendantes, et le résultat distributionnel est approximatif.

Proposition 2.24 (Test de Mardia (1970))

C’est un test basé sur le coefficient d’asymétrie (skewness) et le kurtosis de la loi normale
multivariée. On définit les moments centrés pour X, Y indépendants comme
 3 
> −1
γ1p = E X − µ) Σ (Y − µ )

20
  2 
γ2p = E X − µ ) > Σ −1 ( X − µ ) .

En pratique, on utilise
d2i = ( X i − X̄ )> S−1 ( X i − X̄ )
dij = ( X i − X̄ )> S−1 ( X j − X̄ )
avec les estimateurs
1
b1p =
γ
n2 ∑ ∑ d3ij
i j
1
γ
b2p =
n ∑ d4i
i

Mardia a montré que, pour ν = p( p + 1)( p + 2)/6,

nγ
b1p · 2
∼ χν
6
8p( p + 2)
 
·
b2p ∼ N p( p + 2),
γ
n

Le logiciel R calcule et rapporte les deux statistiques, et indique le rejet dès qu’une des
deux pointe dans ce sens.
Proposition 2.25 (Test de multinormalité basé sur la distance d’énergie)
La distance d’énergie entre deux vecteurs aléatoires X, Y de dimension d est définie par

E ( X, Y ) = 2E| X − Y |d − E| X − X 0 |d − E|Y − Y 0 |d

pour X 0 (Y 0 ) copies iid de X (Y) et | · |d dénote la norme euclidienne. On peut montrer

différentes propriétés liées à cette distance :
1. Elle est invariante par rapport aux rotations
2. E ( X, Y )1/2 définit un métrique sur l’ensemble des lois à d-variables et E ( X, Y ) est
zéro si et seulement si X et Y sont identiquement distribuées.
3. Si d = 1 on a E ( X, Y ) = 2 ( F ( x ) − G ( x ))2 dx où F et G sont les fonctions de
R

répartition de X et Y respectivement.  d +1 
Γ( 2 )
4. Pour une loi Z ∼ Nd (0d , Id ) on a E| Z − Z 0 |d = 2 .
Γ( d2 )

Pour appliquer ces résultats à un test de normalité on commence par standardiser

l’échantillon Y = y via x = S−1/2 (y − ȳ), et on le compare à une variable Z ∼
Nd (0d , Id ). On calcule la statistique
  
d +1
 
2 n Γ 2
  n n
1
n E = n  ∑ E | xi − Z |d − 2     − 2 ∑ ∑ | xi − x k |d 
n i =1 Γ 2 d n i =1 k =1

Le terme E| xi − Y |d est calculable puisque | a − Y |2d suit une loi χ2 non centrée, dont le

21
coefficient de non centralité dépend de a et le nombre de degrés de libertés est aléatoire.
En pratique, on utilise une troncation de la série et la distribution de la statistique est
obtenue par auto-amorçage (bootstrap) paramétrique ou via une approximation asymp-
totique. Le test est consistent sous toutes les alternatives fixes.

2.4.3. Transformations améliorant la normalité

Il existe des techniques basées sur la notion de robustesse. On considère ici une trans-
√
formation, comme par exemple x, log( x ), x −1 , logit( x ), 12 log((1 + x )/(1 − x )), etc. La
dernière est la transformation de Fisher. On a également le cas de la transformation de
Box-Cox, soit
Proposition 2.26 (Transformation de Box-Cox)
On considère la transformation

 X λ −1 si λ 6= 0
X (λ) = λ
log( X ) si λ = 0.

Pour éviter de prendre le logarithme ou une puissance d’un nombre négatif, on peut
considérer X + c. Pour estimer λ, on utilise la log-vraisemblance profilée : pour une
séquence de valeurs λ, on obtient les maximums de vraisemblance de la loi normale
univariée en utilisant comme données xiλ pour i = 1, . . . , n et λ fixe. On peut ainsi
calculer un intervalle de confiance pour λ. Pour tester le paramètre, on peut faire un
test de rapport de vraisemblance,

`0∗ supθ∈Ω0 `
L= ∗ =
`01 supθ∈Ω0 tΩ1 `

·
et −2L ∼ χ2s−r , où s = dim(Ω0 t Ω1 ) et r = dim(Ω0 ).

22
 Chapitre 3
Inférences relatives aux paramètres de
lois normales
3.1 Distribution de X̄ et S

iid
Soit X> = ( X 1 , X 2 , . . . , X n ) où X i ∼ N p (µ, Σ). Alors,

1. la loi de la moyenne arithmétique est X̄ ∼ N p (µ, Σ/n)

2. la loi de la variance empirique est (n − 1)S ∼ W p (n − 1, Σ). On a effectivement

1 >
 
>
( n − 1) S = X I − 11 X
n

où A := I − n1 11> avec rang(A) = n − 1. On a B = n−1 1. Qui est plus, X̄ et S sont

indépendants. On vérifie que A = A2 et que B> A = 0. Ainsi, (n − 1)S ∼ W p (n − 1, Σ)
(et indépendant de X̄).
3. T 2 = n( X̄ − µ)> S−1 ( X̄ − µ) ∼ T p,n
2
−1 . On réécrit cette expression sous la forme

n−1
T2 = n( X̄ − µ)> S−1 ( X̄ − µ)
n−1
√ √ √ √
= (n − 1)( n X̄ − nµ)> ((n − 1)S)−1 ( n X̄ − nµ)
∼ N p (0, Σ)> × W p (n − 1, Σ) × N p (0, Σ)

2 ( n −1) p
suit une loi de Hotelling T p,n −1 ou n− p F p,n− p .

On rappelle le test du rapport des vraisemblances (LRT en anglais) pour un échantillon

X> = ( X 1 , . . . , X n ) dépendant du paramètre θ. On veut tester H0 : θ ∈ Ω0 contre
H1 : θ ∈ Ω1 . On construit λ(x) = L0∗ /L01 ∗ où L∗ = sup
0
∗
θ∈Ω0 L et L01 = supθ∈Ω0 ∪Ω1 L
pour une vraisemblance L.

La région de rejet de H0 (α étant fixé) est définie par R : {x : λ(x) < c} où c est tel que
supθ∈Ω0 Pθ (X ∈ R) = α. On recourt souvent à une approximation.

Théorème 3.1
Soit Ω0 ∪ Ω1 région de Rd avec r = dim(Ω0 ). Alors, sous des conditions suffisantes de
régularité et pour chaque θ ∈ Ω0 , alors la statistique

·
−2 log(λ) ∼ χ2d−r .

autour de b
θ, pour obtenir

23
3.2 Tests relatifs à une moyenne

On s’intéresse à l’hypothèse H0 : µ = µ0 versus H1 : µ 6= µ0 ou plus généralement

Rµ = r versus Rµ 6= r où l’on assume que Σ est connue. Dans ce cas,

n 1   n
` =C− log(|Σ|) − tr Σ−1 S XX − ( X̄ − µ)> Σ−1 ( X̄ − µ)
2 2 2
avec
n 1   n
sup ` = C − log(|Σ|) − tr Σ−1 S XX − ( X̄ − µ0 )> Σ−1 ( X̄ − µ0 )
µ = µ0 2 2 2
n 1  
sup ` = C − log(|Σ|) − tr Σ−1 S XX
µ ∈R p 2 2

pour une constante d’intégration C ∈ R. Le maximum est atteint pour µ

b = X̄ et donc

−2 log(λ) = n( X̄ − µ0 )> Σ−1 ( X̄ − µ0 ) ∼ χ2p

ou plus généralement
  −1
−2 log(λ) = n(R X̄ − r)> RΣR> (R X̄ − r) ∼ χ2q

où q est le rang de r.

On considère maintenant le cas Σ inconnue : la log-vraisemblance est

n n   
` =C− log(|Σ|) − tr Σ−1 n−1 S XX + ( X̄ − µ)( X̄ − µ)> .
2 2
Sous H0 , on a
n   np
sup ` = C − log n−1 S XX + ( X̄ − µ0 )( X̄ − µ0 )> −
µ=µ0 ,Σ 2 2

b 0 = n−1 S XX + ( X̄ − µ )( X̄ − µ )> . Sous l’hypothèse alternative H1 ,

avec Σ 0 0

n np
sup ` = C − log(|S∗ |) −
µ ∈R p 2 2

b = n−1 S XX et µ
avec Σ b = X̄

−2 log(λ) = −n log(|S∗ |) + n log S∗ + ( X̄ − µ0 )( X̄ − µ0 )>

= −n log(|S∗ |) + n log S∗ I + S∗ −1 ( X̄ − µ0 )( X̄ − µ0 )>

 

= n log I + S∗ −1 ( X̄ − µ0 )( X̄ − µ0 )>

= n log 1 + ( X̄ − µ0 )> S∗ −1 ( X̄ − µ0 )
 

en utilisant le résultat A7.10 de l’annexe, qui dit |B + CC> | = |B|(1 + C> B−1 C). On
calcule d = p + p( p + 1)/2 et r = p( p + 1)/2. Cette expression plus petite qu’une

24
constante c ne mène pas au rejet de l’hypothèse nulle H0 si

n log 1 + ( X̄ − µ0 )> S∗ −1 ( X̄ − µ0 ) < c

 

⇔ (n − 1)( X̄ − µ0 )> S∗ −1 ( X̄ − µ0 ) < c0

et donc

n( X̄ − µ0 )> S−1 ( X̄ − µ0 ) ∼ T p,n

2
−1

Pour un test de restriction générale linéaire, on obtient

n(R X̄ − r)> (RSR> )−1 (R X̄ − r) ∼ T 2 (rang(R), n − 1).

Remarque (Intervalles simultanés)

On prend a ∈ R p et on pose Ya = a> X ∼ N (µY , σY2 ) avec Var (Ya ) = a> Σa et E (Ya ) =
a> µ. On utilise la statistique

ȳ − µY a> ( X̄ − µ0 )( X̄ − µ0 )> a
→ t2n = n
a> S∗ a
q
SY2 /n

pour tout a, et on a un quotient de deux formes quadratiques. On peut donc chercher le

maximum maxa t2n et l’on applique le théorème de Cauchy–Schwartz. On rappelle que

(z> y)2 ≤ (z> z)(y> y) (3.1)

avec égalité si et seulement si z ∝ λy. On pose ici z = Σ1/2 a et y = Σ−1/2 ( X̄ − µ). On

obtient

max t2n = n( X̄ − µ0 )> S∗ −1 ( X̄ − µ0 )

a

Si l’on fait l’intervalle de confiance pour toutes les directions possibles et que l’on consi-
dère leur enveloppe, on obtient l’ellipse. Dans la direction X1 , on a, avec probabilité
√
(1 − α), X̄1 ± tn−1,1−α/2 Sn / n . En dimension p, la probabilité de couvrir le rectangle
est (1 − α) p et δ = P (rejeter H0 ) = 1 − (1 − α) p . Ceci illustre la nécessité de faire du
multivariée. La correction de Bonferroni consiste à choisir α tel que δ = 0.05. Ce n’est
pas nécessairement la meilleure région (c’est l’ellipse), mais on obtient un risque global
de 0.05.

Pour l’intervalle de confiance, on a

a> b
max = = b> Σb
a a> Σa

pour a = λΣ−1 b. L’intervalle de Bonferroni est dérivé en utilisant pour m test un seuil

25
de signifiance de αi = α/m, afin que le risque global soit de valeur α au maximum.

Exemple 3.1
On prend a1 = (1, . . . , 0)> , a2 = (0, 1, 0, . . . , 0)> , etc. pour calculer les intervalles de
confiance successivement. On a
s
Si2
X̄i ± t
n n−1,1−α/2

un intervalle de confiance pour µi . La probabilité de couvrir la vraie valeur est 1 − α.

En supposant l’indépendance entre les composantes, on obtient que la probabilité que
les intervalles couvrent les µi est 1 − (1 − α) p . Prenons p = 6 en exemple, on a pour
α = 0.05 que δ = 0.26.

Exemple 3.2
Prenons un exemple numérique concret : on considère un échantillon de n = 10 obser-
vations tirées d’une loi bivariée normale de moyenne arithmétique X̄ = (2.5, 3)> et de
matrice de variance-covariance empirique S = 21 12 . On s’intéresse aux intervalles de

confiance pour µ.

Considérons d’abord l’ellipse correspondant aux intervalles de confiance simultanée.

La matrice de précision est S−1 = −2/3 −1/3
1/3 2/3 . Dans ce cas,

1 18 1 ( n − 1) p
  
( X̄ − µ)> S−1 ( X̄ − µ) ≤ 1.0035 = F2,8 (0.95) = F p,n− p (1 − α).
10 8 n (n − p)

On procède au calcul des valeurs propres et vecteurs propres de S−1 . Le polynôme

caractéristique est (2/3 − λ)2 − 1/9 et en égalant ce dernier à zéro, on obtient λ1 =
√ √
1, λ2 = 1/3. Les vecteurs propres sont donc v1 = (1, −1)> / 2 et v2 = (1, 1)> / 2 et
l’équation de l’ellipse est x12 + x22 /3 = 1.0035. On peut isoler des droites tangentes le
( n −1) p
long de cet ellipse, avec b2 = a> Sa × F
n(n− p) p,n− p
(1 − α ).
r
> > 10 > 10 27 10
a = (1, 1) : a X̄ = , a Sa = 6, a X̄ ± b = ± 4.46 = ± 2.4538
2 2 20 2
r
9
a = (0, 1)> : a> X̄ = 3, a> Sa = 2, a X̄ ± b = 3 ± 4.46 = 3 ± 1.4167
20
r
5 5 9 5
a = (1, 0)> : a> X̄ = , a> Sa = 2, a X̄ ± b = ± 4.46 = ± 1.4167
2 2 20 2
r
1 1 9 1
a = (1, −1)> : a> X̄ = − , a> Sa = 2, a X̄ ± b =− ± 4.46= − ± 1.4167
2 2 20 2

ce qui donne 3.046 ≤ µ1 + µ2 ≤ 7.954, 1.0833 ≤ µ1 ≤ 3.9167, 1.5833 ≤ µ2 ≤ 4.4167 et

finalement −1.9167 ≤ µ1 − µ2 ≤ 0.9167. Si on se restreignait à l’analyse univariée,
on obtiendrait des intervalles de confiance marginaux. Si X ∼ N p (µ, Σ), on a a> X̄ ∼

26
F IGURE 3 – Intervalles de confiance marginaux et simultanés pour le jeu de données
factices. Les intervalles marginaux sont en bleu,les intervalles de Bonferroni en vert et
les projections de l’ellipse de confiance conjointe en pointillés.

5
4
3
µ2

2
1
0

0 1 2 3 4 5

µ1

N a> µ, a> Σa/n et donc

√ a> X̄ − a> µ
n √ ∼ T n −1
a> Sa

et on choisissant a = (1, 0)> ou a = (0, 1)> . On a donc

√ 2.5 − µ
µ1 : 10 √ 1 ≤ 1.833 → 1.68 ≤ µ1 ≤ 3.32
2
√ 3 − µ1
µ2 : 10 √ ≤ 1.833 → 2.18 ≤ µ2 ≤ 3.82,
2

pour un niveau de confiance de 95%, comparativement aux intervalles de Bonferroni,

√
qui sont de la forme X̄i ± Tn−1 (α/2m) Si /n et légèrement plus larges, avec cette fois-
ci 1.3 ≤ µ1 ≤ 3.7 et 1.8 ≤ µ2 ≤ 4.2. Ces deux intervalles sont illustrés par des bandes
bleus et vertes, respectivement, tandis que les axes de l’ellipse sont représentés par
des lignes pleines et les intervalles simultanés à la frontière de l’ellipse par des lignes
pointillées.

27
3.3 Tests relatifs à une variance

Essentiellement, on traitera de 3 hypothèses nulles H0 , à savoir

H0 : Σ = Σ 0 , H0 : Σ = kΣ0 , H0 : Σ12 = 0.

Pour chacune d’elle on part de la log-vraisemblance d’une multinormale, écrite comme

n n   n
` =C− log(|Σ|) − tr Σ−1 S∗ − ( X̄ − µ)> Σ−1 ( X̄ − µ)
2 2 2
On doit maximiser cette quantité sous H0 et sous H0 t H1 . On obtiendra pour chaque
maximisation µ
b = X̄ et donc

n b 0 |) − n tr Σ
 −1 
b 0 S∗
sup ` = C −log(|Σ
H0 2 2
n b 01 |) − n tr Σ
 −1 
b 01 S∗ .
sup ` = C − log(|Σ
H0 t H1 2 2

b 01 = S∗ , on aura
Puisque Σ

−2 log(λ) = npA0 − np log( G0 ) − np,

où A0 est la moyenne arithmétique des valeurs propres de Σ b 0−1 S∗ et G0 la moyenne

géométrique de ces mêmes valeurs propres. On peut maintenant calculer dans chacun
des cas Σb 0.
Exemple 3.3
Sous l’hypothèse nulle H0 : Σ = Σ0 versus H1 : Σ 6= Σ0 on a Σ b 0 = Σ0 ce qui aura
comme conséquence que A0 et G0 sont les moyennes arithmétiques et géométriques de
Σ0−1 S∗ et la statistique suit une loi χ2 avec 21 p( p + 1) degrés de liberté.
Exemple 3.4
k = tr(Σ0−1 S∗ )/p et
Sous l’hypothèse nulle H0 : Σ = kΣ0 versus H1 : Σ 6= Σ0 , on a que b
on obtient finalement

A0
 
·
−2 log(λ) = np log ∼ χ21 p( p+1)−1 .
G0 2

où A0 et G0 sont à nouveau les moyennes arithmétiques et géométriques de Σ0−1 S∗ .

Dans le cas où la matrice Σ0 = I p on fait ce qu’on appelle un test de sphéricité. On peut
aussi généraliser ce test et tester R = I p où R est la matrice des corrélations. Dans ce

1/p
cas, on aura A0 = p/p = 1 puisque diag( R) = 1 p , tandis que G0 = ∏i λi∗ . La
statistique devient alors

·
b |) ∼ χ21
−2 log(λ) = −n log(| R p( p+1)− p
2

28
et suit approximativement une loi khi-deux avec 12 p( p − 1) degrés de liberté. Box (1949)
a démontré que l’approximation est meilleure si on remplace n par n − 1 − 61 (2p + 5).

Exemple 3.5 (Test pour la nullité d’une sous-matrice de variance)

La situation suivante resurgira lors de l’analyse canonique. On veut tester l’hypothèse
H : Σ12 = 0, où µ est inconnue, contre l’alternative Σ12 6= 0. Sous H0 , la vraisemblance
peut se scinder en deux et on aura
! !
∗
S11 0 µ
b1
Σ
b0 =
∗ , b=
µ = X̄.
0 S22 µ
b2

On obtiendra ainsi
! −1
∗ −1 ∗
! !
∗ ∗ ∗
b 0−1 S∗ S11 0 S11 S12 I S11 S12
Σ = ∗ ∗ ∗ = ∗ −1 ∗
0 S22 S21 S22 S22 S21 I

b 0−1 S∗ = p et
avec tr Σ

∗ ∗ ∗ −1 ∗
 −1 
b 0 S∗ = |S∗ | |S22 − S21 S 11 S12 |
det Σ ∗ ∗ = ∗
|S11 ||S22 | |S22 |

et la statistique du rapport de vraisemblance est

|S∗ |
 
−2 log(λ) = n( p1 + p2 ) − n log − n ( p1 + p2 )
| Σ0 |
∗ −1 ∗ ∗ −1 ∗
= −n log I − S22 S21 S11 S12
k
!
= −n log ∏(1 − λ bi)
i =1
·
∼ χ2m

avec
1 1 1
m= ( p1 + p2 )( p1 + p2 + 1) − p1 ( p1 + 1) − p2 ( p2 + 1) = p1 p2
2 2 2
∗ −1 ∗ ∗ −1 ∗ ∗ −1 ∗ ∗ −1 ∗
et 1 − b
λi est une valeur propre de I − S22 S21 S11 S12 ou I − S11 S12 S22 S21 . On
∗ ∗
a cette décomposition en fonction de S11 ou de S22 , ce qui fait que k = min( p1 , p2 ).
Si on travaille plutôt avec la matrice R = D−1/2 S∗ D−1/2 on constate après quelques
calculs que la statistique ne change pas. Idem si l’on utilise l’estimateur sans biais de la
variance, S, plutôt que le maximum de vraisemblance.

On peut montrer que le ratio du déterminant de deux lois Wishart est distribuée exac-
tement comme Λ( p2 , n − 1 − p1 , p1 ), qui est une loi de Wilks Λ( a, b, c). Cette dernière
est liée à la loi F de Fisher pour certaines valeurs des paramètres. Bartlett a proposé

29
une approximation qui est

1
 
− b − ( a − c + 1) log(Λ( a, b, c)) ∼ χ2ac .
2

Ici, ac = p1 p2 = m. Sous H0 : Σ12 = 0, cette approximation donne

 k
1
  
λi ∼ χ2p1 p2 .
− n − ( p1 + p2 + 3) ∑ log 1 − b
2 i =1

Ce test pourrait servir par exemple à vérifier une possible association entre types d’exa-
mens (livre ouvert ou fermé) dans le jeu de données examens.

3.4 Comparaison de deux ou plusieurs populations normales

On considère le test d’égalité des moyennes de populations. Cette situation est une
généralisation de l’ANOVA (analyse de variance), appelée MANOVA.

Dans le cas univarié, on compare les moyennes des deux populations µ1 − µ2 = d0 , à

l’aide du test de Student

X̄ (1) − X̄ (2) − d0
q ∼ Tν
S p n1 + n12
1

où X̄ (k) et Sk sont respectivement la moyenne et la variance du groupe k ∈ {1, 2}, et

(n1 − 1)S12 + (n2 − 1)S22

S2p =
n1 + n2 − 2

est la matrice de variance pour les observations groupées.

Proposition 3.2 (Analyse de variance multivariée : égalité de moyennes)
Soit un échantillon de n observations issues de K populations multinormales. Les don-
nées tirées de la population k seront dénotées X (k) et leurs moyennes et variance em-
piriques par X̄ (k) et Sk . On teste l’hypothèse H0 : µ1 = µ2 = · · · = µK en supposant
Σ1 = Σ2 = · · · = ΣK inconnue, contre l’alternative où il y a au moins une différence
dans les moyennes. On procède au rapport des vraisemblances

1 K    
` = C− ∑ nk log |Σ| + nk tr Σ−1 S∗k + dk d>
k
2 k =1

où dk = X̄ (k) − µk , n = ∑kK=1 nk et

nk
1 (k) (k)
S∗k = ∑ (X j − X̄ (k) )( X j − X̄ (k) )> .
nk j =1

30
 b = S∗ , 3 et
Sous H0 , le maximum de la log-vraisemblance est obtenu pour Σ T

K nk
(k)
bk = X̄ = n−1
µ ∑ ∑ Xj .
k =1 j =1

b = 1 ∑K nk S∗ =
. Sous H0 ∪ H1 , µk = X̄ (k) , et l’estimateur conjoint de la variance est Σ n k =1 k
∗ 4 K ∗ ∗ ∗
S p . On va dénoter Wc=∑ n S
k =1 k k = nS p la variance intra-groupes, T
b = nS T et

K
b=
B ∑ nk (X̄ (k) − X̄ )(X̄ (k) − X̄ )> ,
k =1

la variance inter-groupe. On peut décomposer la variance totale via T b=W c + B.b Cette
dernière est facile à dériver si l’on note qu’en additionnant et en soustrayant X̄ (k) dans
chaque terme,

K nk
(k) (k)
b=
T ∑ ∑ (X j − X̄ )( X j − X̄ )>
k =1 j =1
K nk K nk
(k) (k)
= ∑ ∑ (X j − X̄ (k) )( X j − X̄ (k) )> + ∑ ∑ (X̄ (k) − X̄ )(X̄ (k) − X̄ )>
k =1 j =1 k =1 j =1
K nk K
(k) (k)
= ∑ ∑ (X j − X̄ (k) )( X j − X̄ (k) )> + ∑ nk (X̄ (k) − X̄ )(X̄ (k) − X̄ )>
k =1 j =1 k =1

=W
c + B.
b

Les termes croisés s’annulent puisque

K  nk 
(k) (k)
∑ ∑ j
( X − X̄ ) ( X̄ (k) − X̄ )> = 0.
k =1 j =1

On obtient la statistique
 n
c + n log T


−2 log(λ) = 2 − log W


2 2
b
et donc
W −1
λ2/n = Tb −1 W c −1 B
c
c = = I+W b ∼ Λ( p, n − K, K − 1)
b +W
B c

laquelle loi de Wilks est la distribution exacte si les hypothèses distributionnelles de

départ sont satisfaites (variance conjointe et multinormalité), et suit asymptotiquement

3. Le sous-indice T est utilisé pour dénoter la variance totale, car on peut considérer les observations
comme provenant du même échantillon.
4. L’indice p est en référence au terme anglais pooled, tandis que W dénote within-group variance et B
between-group variance.

31
une loi khi-deux. L’approximation de Bartlett est donnée par

1
 
·
− −n − (K + p + 2) log(Λ) ∼ χ2p(K−1) .
2

Exemple 3.6
Le cas K = 2 permet d’obtenir une expression parlante. On calcule

b = n1 ( X̄ (1) − X̄ )( X̄ (1) − X̄ )> + n2 ( X̄ (2) − X̄ )( X̄ (2) − X̄ )>

B

et X̄ = (n1 X̄ (1) + n2 X̄ (2) )/(n1 + n2 ). En mettant les termes en évidence,

n1 n2
b=
B ( X̄ (1) − X̄ (2) )( X̄ (1) − X̄ (2) )>
n1 + n2

On procède maintenant au calcul de la statistique ; soit

I+W
c −1
b = I + n1 n2 W
B c −1 ( X̄ (1) − X̄ (2) )( X̄ (1) − X̄ (2) )>
n1 + n2
n n c −1 ( X̄ (1) − X̄ (2) )
= 1 + 1 2 ( X̄ (1) − X̄ (2) )> W
n1 + n2

en utilisant un des résultats du résumé d’algèbre linéaire.

Dès lors, avec S p = (n1 S1∗ + n2 S2∗ )/(n1 + n2 − 2) on a

n1 n2 · ( n1 + n2 − 2) p
( X̄ (1) − X̄ (2) )> S− 1
p ( X̄
(1)
− X̄ (2) ) ∼ F
n1 + n2 n1 + n2 − p − 1 p,n1 +n2 − p−1

On rappelle que dans le cas d’une seule population, on avait obtenu en supposant
H0 : µ = µ0 contre l’alternative H1 : µ 6= µ0 , la statistique n( X̄ − µ0 )> S−1 ( X̄ − µ0 ).

3.4.1. Interactions
On considère une ANOVA à deux voies, exprimée comme un modèle linéaire de la
forme

(kg) (kg)
Xi = µ + αk + β g + ε i
(kg) (kg)
Xi = µ + αk + β g + αk : β g + ε i ,

(kg) ind
pour k = 1, . . . , K; g = 1, . . . , G et i = 1, . . . n et ε i ∼ N (0, σ2 ). Ce modèle est appelé
modèle d’analyse de variance (ou ANOVA). Le premier modèle implique que l’effet
quantitatif des facteurs indexés par k et g est purement additif. À l’inverse, si les ef-
fets pour ces facteurs sont plus subtiles, il faudra utiliser un modèle avec interaction,
laquelle est représentée par le terme αk : β g .
Par analogie avec le cas précédent, on a dans le cas multivarié les résultats suivants :

32
Proposition 3.3 (MANOVA à une voie)
Soit le modèle

(k) (k)
Xi = µ + α k + εi

(k) ind
avec εi ∼ N p (0, Σ) pour j = 1, . . . , nk et k = 1, . . . , K. On impose la contrainte
∑kK=1 nk αk = 0. L’hypothèse sous considération est H0 : αk = 0. On peut réécrire

(k) (k)
Xj = X̄ + ( X̄ (k) − X̄ ) + ( X j − X̄ (k) )

et calculer T
b = B c qui donnent λ2/n = |W
b + W, c |/| B c | ∼ Λ( p, n − K, K − 1).
b +W
Asymptotiquement, on aura

K+p
 
− n−1− log(Λ) ∼ χ2p(K −1)
2

sous l’hypothèse nulle.

On a utilisé le critère de Wilks jusqu’à maintenant, mais il y a d’autre critères. Par

exemple, p

• Λ de Wilks : W c / Bb +Wc
• Roy : basé sur la plus grande valeur propre de W c −1 B
b
• Lawley-Hotelling : t = tr(Wc −1 B b)
Pillai : ν = tr B − 1 , réputé moins puissant, mais plus robuste


• b (W + B
c b)

Si k = 2, les quatre critères sont équivalents.

Proposition 3.4 (MANOVA à deux voies avec interactions)

En supposant qu’on a le même nombre d’observations n dans chaque groupe, on pose

(kg) (kg)
Xi = µ + α k + β g + α k : β g + εi

(kg) ind
avec εi ∼ N p (0, Σ) pour k = 1, . . . , K; g = 1, . . . , G et i = 1, . . . n. L’estimateur avec
les données s’exprime sous la forme

(kg) (kg)
Xi = X̄ + ( X̄ (k•) − X̄ ) + ( X̄ (• g) − X̄ ) + ( X̄ (kg) − X̄ (k•) − X̄ (• g) + X̄ ) + ( X i − X̄ (kg) ).

et on décompose la matrice de variance totale en plusieurs composantes,

b=L
T b+C
b + bI + W
c

pour respectivement les effets lignes, colonnes, interactions et la matrice de variance

33
intra-groupe. Les matrices sont définies de façon usuelle, c’est-à-dire

K G n
(kg) (kg)
b=
T ∑ ∑ ∑ (X i − X̄ )( X i − X̄ )>
k =1 g =1 i =1
K
b = nG
L ∑ (X̄ (k•) − X̄ )(X̄ (k•) − X̄ )>
k =1
G
b = nK
C ∑ (X̄ (•g) − X̄ )(X̄ (•g) − X̄ )>
g =1
K G
bI = n ∑ ∑ (X̄ (kg) − X̄ (k•) − X̄ (•g) + X̄ )(X̄ (kg) − X̄ (k•) − X̄ (•g) + X̄ )>
k =1 g =1
K G n
(kg) (kg)
c=
W ∑ ∑ ∑ (X i − X̄ (kg) )( X i − X̄ (kg) )> .
k =1 g =1 i =1

Les hypothèses suivantes peuvent être intéressantes (les résultats sont valides sous
H0 ) :

• H0 : αk : β g = 0 (interactions nulles), calculée avec

W
c
∼ Λ( p, KG (n − 1), (K − 1)( G − 1))
bI + W
c

• H0 : αk = 0 (effets lignes nuls), calculée avec

W
c
∼ Λ( p, KG (n − 1), K − 1)
b +W
L c

• H0 : β g = 0 (effets colonnes nuls), calculée avec

W
c
∼ Λ( p, KG (n − 1), G − 1)
b +W
C c

Si les interactions ne sont pas significatives, et qu’on élimine αk : β g du modèle pour

tester la présence d’effets lignes et colonnes, alors W c On doit toujours inclure
c 7→ bI + W.
les effets principaux pour que le modèle soit sensé.

Proposition 3.5 (Test d’égalité des variances et statistique de Box)

Soit H0 : Σ1 = Σ2 = · · · = ΣK contre l’alternative H1 : Σk 6= Σ j pour au moins un
bk = X̄ (k) , Σ
k 6= j. Sous H0 ∪ H1 , on a µ k bk = X̄ (k) et Σ
b k = S∗ et sous H0 : µ b = S∗p = W
c /n.

34
La statistique est

K |S∗p |
!
·
−2 log(λ) = ∑ nk log ∼ χ2p( p+1)(K−1)/2
k =1
|S∗k |

La statistique de Box (1949) est

K 
|S p |

M=γ ∑ (nk − 1) log |Sk |
k =1
avec
K
2p3 + 3p − 1
!
1 1
γ = 1−
6( p + 1)(K − 1) ∑ nk − 1 − n − K .
k =1

une meilleure approximation à distance finie que la distribution asymptotique.

35
 Chapitre 4
Analyse en composantes principales
L’analyse en composantes principales (ACP) est une technique d’algèbre linéaire qui
remonte à Fisher et Hotelling. La motivation est multiple : d’une part, elle sert à la ré-
duction de la dimension, si par exemple les données vivent dans un sous-espace de
taille plus faible. D’autre part, elle permet d’interpréter les données de façon succincte,
permettant parfois une meilleur perception du problème. Cela peut être utile pour com-
biner des variables mesurées sur des échelles différentes, par exemple les résultats aux
différents épreuves d’un décathlon. Le quotient intellectuel est un exemple d’analyse
en composantes principales (ACP).
Faire une ACP reviendra à diagonaliser la matrice de variance covariance ou la matrice
des corrélations et en extraire vecteurs propres orthonormés et valeurs propres pour
l’analyse.
Définition 4.1 (Analyse en composante principale)
Soit X un vecteur aléatoire de dimension p × 1 avec E ( X ) = µ et Var ( X ) = Σ. La trans-
formation C = P> ( X − µ) où P est une matrice orthogonale p × p telle que P> ΣP = Λ
diagonale (où sans perte de généralité, λ1 ≥ λ2 ≥ . . . ≥ λ p > 0). La ie composante
principale est alors définie par Ci = pi> ( X − µ) où pi est le ie vecteur propre corres-
pondant à la valeur propre λi .

En multipliant une observation centrée par P>, on la transforme en un vecteur de lon-

gueur p dont les composantes sont les composantes principales. En fait, cela revient à
exprimer l’observation dans une autre base.
Proposition 4.2 (Propriétés de l’ACP)
1. E (Ci ) = 0
2. Var (Ci ) = λi et Cov Ci , Cj = 0 pour tout i 6= j. En effet,

 
Var (C ) = Var P> ( X − µ) = P> ΣP = Λ.

p
3. ∑i=1 Var (Ci ) = tr(Σ), appelée inertie totale. On a

p
∑ Var (Ci ) = tr(Λ) = tr(ΛP> P) = tr(PΛP> ) = tr(Σ)
i =1

4. Pour tout a tel que a> a = 1, Var a> X est maximale pour a = p1 et vaut λ1 . En

utilisant les multiplicateurs de Lagrange, on cherche à maximiser Var a> X = a> Σa

et l’on a

    ∂L = 2Σa − 2ξa =0
max L = max a> Σa − ξ a> a − 1 ⇒ ∂a
a a  ∂L = −(a> a − 1) = 0
∂ξ

36
 qui donne Σa = λa et a> a = 1.
5. Pour tout vecteur a tel que a> a = 1 et a> pi = 0 pour i = 1, . . . , j − 1, alors
Var a> X est maximale pour a = p j et vaut λ j .

Remarque
1. Si X ∼ N p (µ, Σ), alors la 1e composante principale s’interprète comme l’axe prin-
cipal de l’ellipsoïde de concentration.
2. La droite passant par µ de direction p1 est la droite de plus grande variance. C’est
aussi la droite qui minimise l’espérance du carré de la distance de X à cette droite.
p
3. ∑ik=1 λi / ∑i=1 λi représente la proportion de l’inertie totale expliquée par les k pre-
mières variables principales.
4. Les composantes principales ne sont pas invariantes par rapport à des change-
ments d’échelle (car il n’existe pas de relations entre les valeurs propres).

Composantes principales et échantillon

Dans le cas d’un échantillon, remplacer µ et Σ par les estimateurs usuels X̄ et S (ou S∗ ).
Si, par exemple, les unités des variables X sont des sec., kg, km et un taux, il est difficile
d’analyser l’unité résultante, puisque

C1 = a1 X1 + a2 X2 + a3 X3 + a4 X4

L’idée dans ces cas est de centrer et réduire (sans unités) et faire l’ACP après, ce qui
revient à faire une ACP sur R au lieu de S.

4.1 Inférence en ACP

Théorème 4.3
Soit X 1 , . . . , X n un échantillon de loi normale N p (µ, Σ) et les estimateurs usuels X̄ , S.
Pour Σ définie positive avec λ1 6= λ2 6= · · · 6= λ p , on a asymptotiquement

1.

2
 
λ ∼ N p λ,
b Λ2 pour λ = ( λ1 , . . . , λ p ) >
n−1

2.

1 λj
 
pi ∼ N p pi , V avec V i = λi ∑ p p> .
n−1 i ( λ j − λ i )2 j j
b
j 6 =i

3. On a également

  ( pi )r , ( p j ) s 1
Cov ( pi )r , ( p j )s = −λi λ j .
( λ i − λ j )2 n − 1

37
 4. {b
λ1 , . . . , b
λ p } est indépendant de { b
p1 , . . . , b
p p }.

On trouvera une démonstration dans un article d’Anderson (1963). Si l’on modifie les
· 2
données (pour améliorer la normalité), on a log(b
λi ) ∼ N log(λi ), n− 1 .

L’analyse en composantes principales sert davantage de technique d’exploration des

données plutôt que de technique inférentielle.

4.1.1. Tests d’hypothèse

On esquisse deux tests utilisés dans le cadre de l’ACP : le premier a pour hypothèse
nulle l’égalité des dernières p − k valeurs propres. Le deuxième test sert à vérifier
si la proportion d’inertie pour les premières valeurs est égale à une constante, soit
p
∑ik=1 λi / ∑i=1 λi = ψ0 donnée.

Proposition 4.4 (Test d’isotropie partielle)

Soit un test d’égalité des p − k dernières valeurs propres, avec

H0 :λk+1 = λk+2 = · · · = λ p
H1 :{∃i, j ∈ {k + 1, . . . , p} : λi 6= λ j , i 6= j}.

Cette hypothèse est utile dans la mesure où si un critère de sélection nous mène à in-
clure la k + 1e composante, alors toutes les autres devraient être incluses dans la mesure
où H0 ne peut être rejetée. Souvent, on débute avec k = 0.

Le test du ratio de vraisemblance est de la forme usuelle, à savoir

−2 log(λ) = np( A0 − log( G0 ) − 1)

où A0 et G0 sont les moyennes arithmétiques et géométriques des valeurs propres de

b 0−1 S∗ et comme à l’accoutumée, Σ
Σ b 0 est l’estimateur de Σ sous H0 . On peut montrer
∗
que, sous H0 , Σb 0 et S ont les mêmes valeurs propres b λi∗ pour i = 1, . . . , k. Donc les
p ∗
valeurs propres de Σ λ∗ , . . . , b
b 0 sont (b
1 λ∗ , b
λ∗ , . . . , b
k λ∗ ) avec b
λ∗ = ∑ λbj /( p − k ). Les
j = k +1
valeurs propres de b 0−1 S∗
Σ sont ainsi

λ∗ λ∗p
!
1, . . . , 1, k+1 , . . . ,
b b
λ∗
b λ∗
b

et donc A0 = 1 et
p−k
1  p
  p−k

(b λ∗k+2 · · · b
λ k +1 b λ∗p ) p−k g0 p
G0 =   = .
b∗
λ λ∗
b

38
Sous H0 , on a k valeurs propres distinctes et 1 valeur propre commune, et k vecteurs
propres avec p composantes, mais k (k + 1)/2 restrictions découlant de l’orthogonalité,
soit k + 1 + kp − k(k + 1)/2 paramètres. Sous H0 t H1 , il y a p( p + 1)/2 contraintes,
d’où

g0
 
·
−2 log(λ) = −n( p − k) log ∼ χ2( p−k+2)( p−k−1)/2 .
λ
b ∗

Proposition 4.5 (Test de proportion d’inertie)

Soit un test d’hypothèse que les p − k composantes représentent moins qu’une certaine
proportion de la variabilité totale,

λ1 + · · · + λ k λ1 + · · · + λ k
H0 : = ψ, H1 : 6= ψ.
λ1 + · · · + λ p λ1 + · · · + λ p

Si b
λi sont les valeurs propres de S, alors ψ λk )/(b
λ1 + · · · + b
b = (b λ p ). En
λ1 + · · · + b
appliquant l’équation A.10.2, on montre que

·
b∼
ψ N (ψ, τ 2 )

avec

2 tr(Σ2 ) λ21 + · · · + λ2k

τ2 = 2
(ψ2 − 2αψ + α2 ), α= .
(n − 1) (tr(Σ)) λ21 + · · · + λ2p

On peut ainsi déterminer un intervalle de confiance asymptotique pour ψ.

Si l’on utilise les variables standardisées et la matrice de corrélation R, on peut tester

l’égalité de toutes les valeurs propres, laquelle égalité tient si et seulement si la matrice
de corrélation est I p . Le test devient

2p + 11
 
·
− n− log(det( R)) ∼ χ2p( p−1)/2 .
6

En revanche, tester que les p − k valeurs propres de la matrice de corrélation sont

égales pose de grosses difficultés. Une proposition de Bartlett est d’utiliser la statis-
tique donnée par −(n − 1)( p − k) log( g0 /λ∗ ), où g0 , λ∗ sont relatifs aux valeurs propres
de R. Cette statistique est à comparer au 1 − α quantile d’une distribution χ2 avec
( p − k + 2)( p − k − 1)/2 degrés de liberté. Notez que cette distribution n’est pas la
distribution asymptotique et mène à un test conservateur.

39
4.2 Interprétation et qualité d’une ACP

1. Corrélations entre les composantes principales et les variables initiales :

p
!
 
Cov Cι , X j = Cov p>
ι X, X j ∑ ( p ι ) k Xk , X j


= Cov
k =1
p
∑ ( pι )k Cov Xk , X j = λ ι ( p ι ) j


=
k =1

puisque Σpι = λι pι et Cov Xk , X j = σkj . De plus, Var (Cι ) = λι et Var X j = σjj , ce

qui donne

λι ( pι ) j
Cor Cι , X j = √ √


λι σjj
√
et si on commence avec R au départ, on a λι ( pι ) j puisque les variance sont égales à 1.
On représente graphiquement cet éléments dans le cercle des corrélations, ici de rayon
unitaire. Si la normalisation des variances n’est pas la même pour chaque variable, on
aura une ellipse puisque la variance de chaque axe sera différente.

F IGURE 4 – Cercle des corrélations

PC2

X4 X1
PC1

X3 X2
X5

La variable X4 est très mal représentée dans le plan des deux composantes principales
puisque leur corrélation est presque nulle. À l’inverse, une portion importante de la
variabilité du point X5 est expliquée par les deux premières composantes.
p
2
Interpréter les oppositions sur chaque axe : on a λι = ∑ j=1 Cor Cι , X j et on appelle
parfois contribution de la variable j à l’axe ι

2
Cor Cι , X j
= ( pι )2j
λι

si l’on travaille à partir de la matrice de corrélations. Ce cercle des corrélations est le

pendant pour les variables de la projection des variables. Dans le biplot, les deux gra-
phiques sont superposés. Attention : la représentation doit s’utiliser avec précaution.
On ne peut pas interpréter la proximité des points (individus), mais leur direction

40
puisque souvent on a pas les mêmes unités et les échelles des deux représentations
sont différentes.
2. Effet de “taille” : Si σij > 0 pour tout i et j, on peut montrer qu’alors tous les
coefficients du premier vecteur propre sont de même signe (théorème de Frobenius).
C1 = a1> X est une sorte de moyenne pondérée (toutes les variables n’ont pas le même
poids). Ainsi, la 1e composante principale mesure souvent un effet de “taille”. On
peut aussi remarque que la 2e composante principale, C2 , représente fréquemment un
effet de “forme”.
3. Importance des individus : dans le cas d’un échantillon X de taille n × p, on a C1i =
p1> ( X i• − X̄ ), la projection du ie individu sur la première composante principale. La
variance des points projetés est n1 ∑in=1 Cιi2 = b λι et n−1 Cιi2 /b
λι est la contribution de
l’individu i à la composante ι.
Théorème 4.6 (Frobenius)
Soit A une matrice symétrique avec a jk > 0, alors les composantes du premier vecteur
propre de A sont de même signe.

4. Variables et individus supplémentaires : On essaie parfois de garder des variables

et/ou des individus en réserve pour valider ou confirmer son interprétation. Ceci
peut être une aide.
5. Qualité des représentations :
• Mesure globale : pourcentage d’inertie expliquée par les composantes retenues
p
• Qualité de la représentation des individus : COι = Cιi2 / ∑k=1 Cki 2 et QLT = CO +
1
CO2 . L’indicateur QLT est la qualité de la représentation d l’individu i sur les deux
premières composantes.
6. Nombre d’axes (facteurs, composantes) à retenir :
(a) Critère de Kaiser : en travaillant sur R, ne retenir que les composantes avec valeurs
propres plus grandes que 1.
(b) Critère de Jolliffe : prendre les composantes jusqu’à obtenir plus de 80% de l’iner-
tie totale
(c) Critère du coude (critère graphique) : on trace l’index des valeurs propres contre
λi
(d) Scree test de Cattell : calculer λ1 − λ2 = ε 1 , λ2 − λ3 = ε 2 , . . . et calculer par la
suite ε 1 − ε 2 = δ1 , ε 2 − ε 3 = δ2 . Ne retenir que les valeurs propres λ1 , . . . , λk+1 avec
δ1 , . . . , δk ≥ 0.

41
4.3 ACP, décomposition en valeurs singulières et régression

Soit X la matrice n × p formée des observations X = ( X 1> · · · X > > p

n ) , avec X i ∈ R
centrées (donc assumant que la moyenne des colonnes est égale à zéro). On considère
la sélection de k < p combinaisons linéaires XA ∈ Rk qui, dans un sens, représentent
le mieux les données originales.
La décomposition en valeurs singulière de X est

X = U DV >

où
• U est une matrice orthonormale n × p ; les colonnes de U engendre l’espace des
colonnes de X
• V est une matrice orthogonale p × p ; les colonnes de V engendre l’espace des lignes
de X
• D est une matrice diagonale p × p dont les éléments di sont les valeurs singulières
de X, avec d1 ≥ d2 ≥ · · · d p ≥ 0 et rang(X) = rang( D ).
Une autre formulation de la décomposition singulière prend D comme matrice n × p
contenant le carré des valeurs singulières, à savoir les valeurs propres de XX> et de
X> X, U la matrice n × n des vecteurs singuliers gauches, égale à la matrice de vecteurs
propres de XX> et V > matrice p × p de vecteurs singuliers droits, soit les vecteurs
propres de X> X.
Les composantes principales sont alors les colonnes de XV . On a

(n − 1)S = X> X = V D > U > U DV > = V D2 V >

puisque U est orthonormale et D est diagonale. Dans l’ACP, S = PΛP> . On a donc

P = V et Λ = (n − 1)−1 D2 .

Régression

On rappelle quelques notions de régression linéaire. Le but de cette section est d’uti-
liser la décomposition en valeurs singulières (et l’ACP) pour mieux comprendre les
estimateurs à rétrécissement. On conclura avec deux méthodes qui utilisent une trans-
formation de la matrice de régresseurs X pour la régression.
Dans la régression par moindre carrés, on postule un modèle linéaire de la forme Y =
Xβ + ε, où ε est une matrice d’aléas indépendants et de moyenne nulle. L’estimateur β
b
qui minimise l’erreur moyenne quadratique est la solution du problème

−1 >
min(y − Xβ)> (y − Xβ) ⇒ β
b >
MC = ( X X ) X y si Var (ε) ∝ In
β

min(y − Xβ)> Σ−1 (y − Xβ) ⇒ β

b > −1 −1 > −1
MCP = ( X Σ X ) X Σ y si Var (ε) = σIn
β

42
pour σ connu. Le dernier cas (données hétéroscédastiques) est appelé moindres carrés
pondérés.

Nous considérerons subséquemment les seuls moindres carrés ordinaires, dont l’estimé
sera dénoté β.
b Les valeurs ajustées sont y b = Xβ b = HX y = X(X> X)−1 X> y et par
construction, est orthogonal à l’estimé des résidus b
ε = (I − HX )y.

Si X = U DV> , on obtient

b = U DV > (V D > U > U DV > )−1 V D > U > y = UU > y

b = Xβ
y

où U > y sont les coordonnées de y dans la base orthonormale U. Évidemment, le ré-

sultat découle du fait que U > U = V > V = V V > = I p , D = D > et (V D2 V > )−1 =
V D −2 V > . Puisque X β
b = XV V > β,
b les coefficients par rapport à la matrice de régres-
seurs XV sont orthogonaux et peuvent être obtenus par des régressions simples suc-
cessives. L’exclusion de colonnes de XV n’affectent pas les estimés.

Régression avec pénalité et estimateurs à rétrécissement

Proposition 4.7 (Régression ridge)
Les méthodes de rétrécissement introduisent du biais dans l’estimateur de β au profit
d’une réduction de la variance. Certaines pénalités permettent même de faire de la
sélection de variables. La régression de crête (ridge) ajoute une pénalité sur la norme l2
du vecteur de coefficient, via

min(y − Xβ)> (y − Xβ) + θβ> β ⇒ b = (X> X + θI p )−1 X> y

β r
β

si les erreurs sont homoscédastiques et indépendantes. L’estimateur ridge s’exprime

dans la base de vecteurs propres via

b = (X> X + θI p )−1 X> y

β r
= (V D2 V > + θI p )−1 V DU > y
= V ( D2 + θI p )−1 DU > y
p
dj
= ∑ v j d2 + θ u>j y,
j =1 j

d’où
p d2j
b = Xβ
y b =
r ∑ u j d2 + θ u>j y.
j =1 j

Le coefficient d2j /(d2j + θ ) est un coefficient de rétrécissement.

Proposition 4.8 (Régression LASSO)

Une des pénalités les plus connues est celle de la régression lasso, qui donne le problème

43
d’optimisation quadratique

p
min(y − Xβ)> (y − Xβ) + θ
β
∑ | β j |.
j =1

La solution est obtenable par le biais d’algorithmes dédiés qui permettent de trouver
l’ensemble des solutions en fonction de θ, et ce au même coût de calcul que pour la
régression ridge. La nature de la pénalité pousse les estimés pour certains régresseurs
à zéro, ce qui crée un modèle plus parcimonieux. Là où la pénalité de crête rétrécie de
façon proportionnelle les estimateurs de moindre carré, βbrj = βbj /(1 + θ ), la régression
lasso donne βblasso = sign( βbj ) min{| βbj | − θ, 0}. Mentionnons finalement la pénalité du
j
p
filet élastique, θ ∑ j=1 (αβ2j + (1 − α)| β j |, un compromis entre lasso et ridge.

Les pénalités présentées jusqu’à maintenant ont un désavantage : elles atténuent le

signal même si les coefficients sont larges.
Proposition 4.9 (Régression SCAD)
La pénalité SCAD, de paramètres a > 2 et θ > 0 est une spline quadratique avec noeuds
à θ et aθ de la forme

θ|β | si | β j | ≤ θ
 2j



β j −2aθ | β|+θ 2
p( a, θ ) = − si θ < | β j | ≤ aθ
 2( a −1)

( a + 1 ) θ 2
si | β | > aθ,

j

2

qui donne comme solution si les régresseurs sont orthogonaux


sign( β j ) min{| β j | − θ, 0} si | βbj | < 2θ

 b b
βbscad
j = {( a − 1) βb) j − sign( βbj )) aθ }/( a − 2) si 2θ < | βbj | ≤ aθ

si | βbj | > aθ.

βj
b

Un algorithme de descente par coordonnée est nécessaire pour obtenir une solution.
Comme pour le lasso et la régression ridge, les paramètres optimaux pour la pénalité
sont souvent dérivés par validation croisée. En revanche, les erreurs standards rap-
portées par le progiciel ne prennent pas en compte la sélection de variable et les tests
d’hypothèse basés sur ces dernières ne sont pas valides.

Régression avec dérivés de la matrice de régresseurs

Soit X un tableau n × p non-stochastique centré et réduit (dont la moyenne des colonnes

est de 0 et leur variance de 1) et un vecteur réponse y tel que ȳ = 0, avec le modèle

Y = Xβ + ε, ε ∼ N n (0n , σ 2 I n ).

44
On considère un changement de base de X en prenant zi = X pi , où i = 1, . . . , k ≤ p
Une régression par moindre carrés minimise la distance euclidienne verticale entre les
observations et la droite de régression, tandis qu’une régression avec les composantes
principales minimisera la distance perpendiculaire à la droite.
Proposition 4.10 (Régression sur composantes principales)
On fait la régression de y sur les vecteurs colonnes z1 , . . . , zk , lesquels sont orthogonaux
par construction. Les estimés sont donnés par

k
brcp = ȳ1 + ∑ γ
y bi zi , bi = (zi> zi )−1 zi> y.
γ
i =1

On appelle matrice de chargement la matrice des vecteurs propres de Σ, divisés par

la racine des valeurs propres respectives. Si Σ = PΛP> , la matrice de chargement
est L = PΛ1/2 . Si L1:k = (l1 , . . . , lk ) dénote la sous-matrice de rang k, la meilleure
approximation de la matrice de covariance (ou de corrélation) est donnée par L1:k L1:k > .

Dans le cadre de la décomposition spectrale, U D est une matrice de scores pour la

régression. 5

Le rétrécissement se fait à l’aide du choix de k. En pratique, il est peu courant d’utiliser

la régression en composantes principales puisque cette dernière tente de représenter
de façon parsimonieuse la variation en X sans considération de son pouvoir explicatif.
Toutes les variables de X sont également nécessaires pour les prédictions. Pour illustrer
cette réalité, supposons un instant que la variance σ2 , les coefficients γ et les valeurs
propres λ sont connus.
Proposition 4.11 (Choix optimal de k pour RCP)
Pour minimiser l’erreur moyenne quadratique, on conserve une composante principale
si le carré du biais est plus important que la réduction de variance. L’omission d’une
composante i introduit un biais de γi (puisque l’estimateur γbi est sans biais pour γi et
que les régresseurs zi sont orthogonaux). La réduction de la variance due à l’omission
d’une composante est σ2 ( Z > Z )ii−1 = σ2 /((n − 1)λi ).

Les variables à conserver ne sont donc pas nécessairement les premières composantes
et le choix optimal pour k dépend des paramètres γ, lesquels sont inconnus en pratique.

Régression par moindres carrés partiels (PLS)

Considérons y ∈ Rn et X ∈ Rn× p des tableaux centrés (chaque colonne de X de même

que y ont moyenne zéro). La régression par moindre carrés partiels est une technique
analogue à la régression en composantes principales qui utilise des combinaisons li-
néaires de X et permet de traiter le cas où p  n. Cette technique est particulièrement
5. Notez au passage que R utilise le terme loadings pour dénoter les axes principaux (un abus de lan-
gage).

45
utilisée en chimiométrie. La régression des moindres carrés partiels construit itérati-
vement une base de k ≤ min(n − 1, p) vecteurs de poids Wk = (w1 , . . . , wk ) pour
les colonnes de X, produisant une matrice de scores n × k donnée par l’expression
T = XW.
La première direction est obtenue en maximisant la covariance Cov ( Xw, y), ou bien
encore Cor ( Xw, y) Var ( Xw) pour w ∈ R p sujet à kwk = 1. Ce problème d’optimi-
p

sation contrainte peut être résolu à l’aide des multiplicateurs de Lagrange et on trouve
que w1 est le vecteur propre associé à la plus grande valeur propre de X > yy> X, soit
w1 = X > y/k X > yk. Pour les directions subséquentes, le coefficient wi est la solution
du problème maxw Cov(y, Xw) sujet aux contraintes kwk = 1 et w> Swl = 0 pour
l = 1, . . . , i. La régression PLS effectue un compromis entre l’analyse en composantes
principales (qui maximise la variance) et la régression multiple (l’analyse canonique
dans le cas où Y est n × m), qui maximise la corrélation.
Mathématiquement, les moindres carrés partiels reviennent à minimiser

b = arg min (y − X β)> (y − X β)

β i
β∈Ki (A, b)

avec A = n−1 X > X, b = n−1 X > y, où Ki (A, b) est un espace de Krylov, c’est-à-dire
l’espace engendré par vect{A j−1 b}ij=1 .
Pour la résolution, on exprime le système de régression

X = TP> + E
Y = UQ> + F

où T, U sont les matrices n × k de score, P et Q sont des matrices de chargements de

taille respective p × k et 1 × k, E et F des matrice de résidus, c> ∈ Rk est un vecteur de
contributions et finalement W est une matrice de poids de taille n × k. 6
La résolution du problème des moindres carrés partiels passe par des algorithmes ité-
ratifs, donc l’algorithme NIPALS. Ce dernier alterne entre régression et projection et on
calcule itérativement les colonnes du système. Par exemple, la matrice Tl à l’étape l est
formée par la concaténation des colonnes t1 , . . . , tl , lesquels t sont calculés récursive-
ment. Des généralisation multivariées existent et diffèrent légèrement de la présenta-
tion ci-dessous.
Algorithme 4.1 (Moindres carrés partiels itératif non-linéaire (NIPALS) — univarié)
Soit X et y centrés et le nombre de composantes PLS k :
Initialisation: E0 = X et u0 = F0 = y
1: pour i ← 1, . . . , k

6. Dans le cas multivarié avec Y ∈ Rn×m , on cherche à maximiser Cov ( Xw, Yc) avec comme contrainte
additionnelle kck = 1, etc. C’est donc désormais une matrice de vecteurs orthogonaux C k × m que l’on
cherche.

46
 2: wi ← Ei>−1 Fi−1 /(Fi>−1 Fi−1 ) . régression
3: wi ← wi / k wi k
4: ti ← Ei−1 wi /wi> wi
5: pi ← Ei>−1 ti /ti> ti . régression
6: Ei ← Ei−1 − ti pi> . projection
7: ci ← Fi>−1 ti /ti> ti . régression
8: ui ← Fi−1 /ci
9: Fi ← Fi−1 − ti ci> . projection

Dans le cas de données manquantes, les quantités sont calculées avec les données com-
plètes. Si aucune donnée n’est manquante, on peut remplacer les deux premières étapes
(lignes 2-3) par wi ← Ei>−1 Fi−1 /kEi>−1 Fi−1 k et 4 par ti ← Ei−1 wi .

Les valeurs ajustées après k itérations sont calculées en utilisant la récursion

k ui> ti
b(k) = ȳ + ∑ ti b
y αi , αi = .
ti> ti
b
i =1

On peut réexprimer les coefficients de la régression en fonction de la matrice X de

départ pour le modèle de régression y = XB + F. On a alors

−1
B(k) = Wk ( P>
k Wk ) ck .

Si k = p, l’estimé B(k) est égal à celui de la régression des moindres carrés, tout comme
dans la régression en composantes principales. Si la matrice de régresseurs X est or-
thogonale, les moindres carrés partiels retournent les estimateurs des moindres carrés
après la première itération ; aucune étape subséquente ne modifie les estimés initiaux.
Les erreurs à l’étape h, bfh sont recalculées comme bfh = y − yb(h) .

On peut calculer, comme dans le cas de l’ACP, la proportion de variance expliquée pour
y et X avec h chargements :


Var y
b(h) SCh (y)
pVarh (y) = =
Var (y) y> y
t> >
h th ph ph SCh ( X )
pVarh ( X ) = = .
tr( X > X ) tr( X > X )

Puisque les scores sont orthogonaux, la proportion cumulative des k premières va-
riables PLS est la somme des contributions ∑kh=1 pVarh .

Il peut être difficile de choisir le nombre de composantes pour la régression PLS. Les
principaux critères proviennent de la validation croisée : on peut utiliser la méthode du
canif (jacknife) en laissant de côté une observation et en estimant le modèle sans cette

47
dernière. La somme des carrés des résidus non-standardisée est comme à l’accoutumée
n 2
>

RSSh = bf h bf h = ∑ yi − yb(h)i
i =1
tandis que la somme du carré des résidus prédits par le modèle estimé, dénotée PRESSh
est
n  2
i
PRESSh = ∑ yi − yb− ( h )i
i =1

La statistique PRESSk n’existe pas nécessairement (si par exemple k = n). On définit
quelques statistiques

PRESSh
Q2 ( h ) = 1 −
RSSh−1
h
PRESSh
Q2cum (h) = 1 − ∏
i =1
RSSh−1
PRESSh
R2VC (h) = 1 −
RSS0

avec RSS0 = ∑in=1 (yi − ȳ)2 . En ce qui a trait aux critères de décisions, on a plusieurs
choix

• Ajouter le chargement h si Q2cum (h) > 0.5Q2cum (h − 1)

√
• Ajouter le chargement h si PRESSh ≤ γ RSSh−1 avec γ = 0.95 pour n < 100 et
p

γ = 1 si n ≥ 100.

48
 Chapitre 5
Analyse canonique
Dans un modèle de régression multiple, on pose comme modèle pour une observa-
tion Y = a> X tel que Cor Y, a> X est maximale. Dans le cas où Y est un vecteur

réponse, l’analyse canonique consiste à poser b> Y expliquée par a> X de telle sorte
que Cor b> Y, a> X soit maximale.

5.1 Définition et propriétés de l’analyse canonique

Soit ( X, Y )> un vecteur de composantes X de dimension p1 et Y de dimension p2 ,

p := > :
avec  p1 + p2 , de moyenne (µ1 , µ2 ) = µ et matrice de variance covariance Σ =
Σ11 Σ12

Σ21 Σ22 . On a

 
Cov a> X, b> Y = a> Σ12 b
 
Var a> X = a> Σ11 a
 
Var b> Y = b> Σ22 b.
Ainsi,
  a> Σ12 b
Cor a> X, b> Y = p (5.2)
a> Σ11 ab> Σ22 b

Maximiser cette corrélation donnera une valeur positive (pourquoi ?). 7 De plus, il y
aura un problème de normalisation de a et b. On imposera des contraintes pour maxi-
miser équation (5.2), soit a> Σ11 a = 1 et b> Σ22 b = 1. On résoudra le problème d’opti-
misation à l’aide des multiplicateurs de Lagrange

1 1
L = a> Σ12 b − δ(a> Σ11 a − 1) − γ(b> Σ22 b − 1).
2 2
En dérivant, on obtient

∂L
= Σ12 b − δΣ11 a = 0 ⇒ Σ12 b = δΣ11 a (5.3)
∂a
∂L
= Σ21 a − δΣ22 b = 0 ⇒ Σ21 a = γΣ22 b. (5.4)
∂b
De la première expression, on obtient en inversant les sous-matrices de covariance
−1
Σ11 Σ12 b = δa des expressions pour M1 et M2 , définies comme

−1 −1
M1 a := Σ11 Σ12 Σ22 Σ21 a = δγa = λa
−1 −1
M2 b:= Σ22 Σ21 Σ11 Σ12 b= δγb = λb

7. Si on trouve a telle que la corrélation est négative, prendre −a donnera une valeur positive.

49
Que vaut Cor a> X, b> Y au maximum ? En multipliant par a> et par b> respective-

ment les équations (5.3) et (5.4), on obtient δ = γ = a> Σ12 b. On a donc

a> Σ12 b p
max p = λmax .
a,b a> Σ11 ab> Σ22 b
√ √
La valeur maximale de la corrélation est donc λmax où λmax est la plus grande
valeur propre de M1 et M2 . Les valeurs propres non nulles sont en quantité égales,
soit k = min{ p1 , p2 }. Noter que l’on a rencontré ces matrices lors du test d’hypothèse
H0 : Σ12 = 0. Le problème dans le développement ci-dessus est que ni M1 , ni M2
ne sont symétriques. Il est difficile dans ce cas de trouver les valeurs propres et les
vecteurs propres. En pratique, il n’y a pas de différence que l’on travaille sur R ou
S ; en écrivant R = D −1/2 SD −1/2 pour D la matrice diagonale des variances que l’on
sépare, on obtient les mêmes valeurs propres. On résout le problème en introduisant
les matrices

−1/2 −1 −1/2
N1 = Σ11 Σ12 Σ22 Σ21 Σ11 = KK>
−1/2 −1 −1/2
N2 = Σ22 Σ21 Σ11 Σ12 Σ22 = K> K

−1/2 −1/2
où K = Σ11 Σ12 Σ22 . On montre que M1 , M2 , N1 , N2 ont les mêmes valeurs propres
non-nulles λ1 > λ2 > · · · > λk . On dénote par αi (respectivement βi ) les vecteurs
−1 −1
propres de N1 (respectivement N2 ). Soit a vecteur propre de Σ11 Σ12 Σ22 Σ21 et b vecteur
−1 −1
propre de Σ22 Σ21 Σ11 Σ12 = M2 . Alors, on peut lier les deux par le biais des relations
−1/2
αi> αi = 1, αi> α j = 0, ai = Σ11 αi et ai> Σ11 ai = 1. Pareillement, on a β> >
j β j = 1, βi β j =
−1/2
√
0, b j = Σ22 β j et b>
j Σ22 b j = 1. Notez que k ai k 6 = 1. λi est appelé le ie coefficient de
corrélation canonique, ai> X et bi> Y sont appelés les ie variables canoniques.

Proposition 5.1 (Propriétés)

Posons ηi = ai> X et φi = bi Y ; alors,

1. E (ηi ) = ai> µ1 et E (φi ) = bi> µ2

2. Var (ηi ) = ai> Σ11 ai = αi> αi = 1 et Var (φi ) = bi> Σ22 bi = βi> βi = 1
√
3. Cor ηi , η j = Cor φi , φj = 0 si i 6= j et Cor (ηi , φi ) = ai> Σ12 bi = λi . Ceci repré-

sente la corrélation de la ie variable canonique en X et de la ie variable canonique en Y.

On a également

−1/2 −1/2
Cor ηi , φj = Cov ηi , φj = ai> Σ12 b j = αi > Σ11 Σ12 Σ22 β j = αi > Kβ j

On sait que β j (respectivement αi ) est vecteur propre de K> K (KK> ). On pose donc
K> Kβ j = λ j β j , ce qui implique que KK> Kβ j = KK> dα j = λ j dα j . En notant que
−2 > > −2 >
√
α>
j α j = d β j K K = d λ j β j β j , on vérifie que d = λ j et il en découle que
>
√
Cor ηi , φj = αi λ j α j = 0.

50
Analyse canonique sur échantillon

On remplace Σij par Sij ou Sij∗ . La division par n − 1 ou n est sans importance puisque
−1 −1
ces constantes s’annulent dans le produit Σ11 Σ12 Σ22 Σ21 . Posons R = D −1/2 SD −1/2 ;
on a
−1  −1 
S21 = D11/2 R11 D1/2

D1/2

D1/2 1/2
D1/2 1/2

−1 −1 1/2
S11 S12 S22 1 1 R12 D 2 R22 D2 2 R21 D1

= D1−1/2 R11
−1 −1
R12 R22 R21 D1/2
1

et donc

D1−1/2 R11
−1 −1
R12 R22 R21 D1/2
1 v = λv
−1
R11 −1
R12 R22 R21 D1/2 1/2
1 v = λD1 v.

Notez que la valeur propre maximale est la même si on utilise R ou S. Ainsi, le résultat
de l’ analyse canonique sur R ou sur S est identique.

5.2 Inférence en analyse canonique

On présente quelques tests d’hypothèse associés à l’analyse canonique.

(a) H0 : Σ12 = 0 contre H1 : Σ12 6= 0. On a traité ce test précédemment ; dans le présent

contexte, si S12 ≈ 0, ça ne vaut pas la peine de faire une analyse canonique. La ratio
des log-vraisemblances est

k
!
−2 log(λ) = −n log ∏ (1 − λ i )
i =1

−1 ·
avec λi valeur propre de S22 S21 S11 S12 et la statistique de Wilks donnée par λ2/n ∼
Λ( p2 , n − 1 − p1 , p1 ). On peut utiliser l’approximation de Bartlett, soit
!H
k 0
1
 
·
V = − n − ( p1 + p2 + 3) log ∏ (1 − λ i ) ∼ χ2p1 p2
2 i =1

(b) H0 coefficients de corrélation canonique sont non nuls (ri = 0), où ri2 = λi . La
statistique

1
   
Vj = − n − ( p1 + p2 + 3) log 1 − r2j
2

Le nombre de degrés de liberté est ( p1 − 1)( p2 − 1) pour la statistique V − V1 pour

tester l’hypothèse voulant que les variables canoniques 2 à k ne soient pas significa-
tives, ( p1 − 2)( p2 − 2) pour la statistique V − V1 − V2 pour tester l’hypothèse nulle
que les variables canoniques 3 à k ne sont pas significatives, etc.

51
Interprétation et généralisation
Soit les matrices X et Y, assumées centrées. Les vecteurs a j et b j correspondent aux je
vecteurs de corrélation canonique et Xa j et Yb j représentent les scores des n individus
sur la je variable canonique, tandis que η j = a> >
j x et φ j = b j y représentent la valeur
particulière des scores pour ( x, p
y). Côté représentation graphique de (η1 , φ1 ), on a la
droite correspondante de pente b λi , soit
q
bi =
φ λ i ηi .
b

Si la corrélation canonique est près de 1, alors la relation est presque linéaire de pente
1. On peut tracer également (ηb1 , ηb2 ) versus (φ
b1 , φ
b2 ); il y a forte ressemblance entre les
nuages de points si les corrélations canoniques sont fortes.

5.3 Généralisation : analyse canonique et tableaux de contingence

Tableau 3 – Exemple de tableau de contingence

Cheveux
clairs foncés
bleus 4 0
Yeux verts 1 2
bruns 0 3

On fabrique un tableau X et Y à partir du tableau de contingence créé précédemment,

soit
   
1 0 0 1 0
1 0 0 1 0
   
1 0 0 1 0
   
   
1 0 0 1 0
     
4 0
0 1 0 1 0
   
X= 0 1 0 , Y= 0 1 , X > Y = 1 2
   

0 3
   
0 1 0 0 1
   
0 0 1 0 1
   
   
0 0 1 0 1
   
0 0 1 0 1

Pour faire une analyse canonique, on doit calculer la matrice S = SS11 SS22 12
. En re-


21
vanche, il y a des dépendances ici : la somme des colonnes de X et Y égale 1n et de ce
fait S11 et S22 calculées en recentrant les données ne seront pas inversibles. Les solutions
sont les suivantes
1. Enlever une colonne à X et à Y et effectuer l’analyse canonique sur le résultat par

52
 la suite.
2. Effectuer l’analyse canonique sur
!
e= 1 X >X X >Y
S .
n Y >X Y >Y

Cela donnera une valeur propre nulle, à ne pas considérer. Cela revient (en ne centrant
pas) à trouver une inverse généralisé et conduit à la même solution pour les valeurs
propres non triviales.

Ce résultat correspond exactement à la technique que l’école française de la statistique

appelle l’analyse des correspondances. Elle est due à un algébriste du nom de Benzécri
et sert principalement lors de questionnaires.

Exemple 5.1 (Origine sociale des étudiants)

Une enquête de l’INSEE recense la branche d’étude de 10000 étudiants en relation avec
la catégorie socioprofessionnelle du père. Le résultat de l’enquête est comptabilisé dans
le tableau de contingence.

Tableau 4 – Origine socioprofessionnelle et domaine d’étude d’étudiants, 1975-1976

DROI SEC LETT SCIE MEDD PHAR PLUR IUT Total

EXPA 80 36 134 99 65 28 11 58 511
SALA 6 2 15 6 4 1 1 4 39
PATR 168 74 312 137 208 53 21 62 1035
LICS 470 191 806 400 876 164 45 79 3031
CADM 236 99 493 264 281 56 36 87 1552
EMPL 145 52 281 133 135 30 20 54 850
OUVR 166 64 401 193 127 23 28 129 1131
SERV 16 6 27 11 8 2 2 8 80
AUTR 305 115 624 247 301 47 42 90 1771
Total 1592 639 3093 1490 2005 404 206 571 10000

Les colonnes sont respectivement DROI (droit), SEC (sciences économiques), LETT (Lettres), SCIE (sciences),
MEDD (médecine et dentaire), PHAR (pharmacie) et PLUR (pluridisciplinaire). Les catégories pour les lignes
sont respectivement EXPA (exploitant agricole), SALA (salarié agricole), PATR (patron), LICS (profession
libérale, cadre supérieur), CADM (cadre moyen), EMPL (employé), OUVR (ouvrier), SERV (personnel de
service) et AUTR (autres).

La question d’intérêt est : quelles études poursuivent les enfants des diverses catégories
socio-professionnelles ?

On peut analyser les proportions marginales

Le tableau des profils associés aux lignes (colonnes) permet de voir les proportions.

La moyenne des profils-lignes (l’effectif pondéré par la somme des effectifs marginaux

53
 Tableau 5 – Proportion marginale associée

(a) Lignes (b) Colonnes

EXPA 0.051 DROI 0.16

SALA 0.004 SEC 0.06
PATR 0.103 LETT 0.31
LICS 0.303 SCIE 0.15
CADM 0.155 MEDD 0.20
EMPL 0.085 PHAR 0.04
OUVR 0.113 PLUR 0.02
SERV 0.008 IUT 0.06
AUTR 0.177

Tableau 6 – Profil associés aux lignes

DROI SEC LETT SCIE MEDD PHAR PLUR IUT

EXPA 0.1566 0.0705 0.2622 0.1937 0.1272 0.0548 0.0215 0.1135
SALA 0.1538 0.0513 0.3846 0.1538 0.1026 0.0256 0.0256 0.1026
PATR 0.1623 0.0715 0.3014 0.1324 0.2010 0.0512 0.0203 0.0599
LICS 0.1551 0.0630 0.2659 0.1320 0.2890 0.0541 0.0148 0.0261
CADM 0.1521 0.0638 0.3177 0.1701 0.1811 0.0361 0.0232 0.0561
EMPL 0.1706 0.0612 0.3306 0.1565 0.1588 0.0353 0.0235 0.0635
OUVR 0.1468 0.0566 0.3546 0.1706 0.1123 0.0203 0.0248 0.1141
SERV 0.2000 0.0750 0.3375 0.1375 0.1000 0.0250 0.0250 0.1000
AUTR 0.1722 0.0649 0.3523 0.1395 0.1700 0.0265 0.0237 0.0508

des lignes), est le profil marginal des colonnes :

r nij  ni•  n• j
∑ ni • n
=
n
i =1

On peut diagonaliser le produit de ces deux tableaux de profils ; la décomposition spec-

trale donne les valeurs propres et le pourcentage d’inertie associé.

On choisit deux composantes pour le reste de l’analyse.

La représentation simultanée est alors la suivante :

Comme dans l’analyse canonique et l’analyse en composantes principales, il faut porter

Tableau 7 – Profil associés aux colonnes

EXPA SALA PATR LICS CADM EMPL OUVR SERV AUTR

DROI 0.0503 0.0038 0.1055 0.2952 0.1482 0.0911 0.1043 0.0101 0.1916
SEC 0.0563 0.0031 0.1158 0.2989 0.1549 0.0814 0.1002 0.0094 0.1800
LETT 0.0433 0.0048 0.1009 0.2606 0.1594 0.0909 0.1296 0.0087 0.2017
SCIE 0.0664 0.0040 0.0919 0.2685 0.1772 0.0893 0.1295 0.0074 0.1658
MEDD 0.0324 0.0020 0.1037 0.4369 0.1401 0.0673 0.0633 0.0040 0.1501
PHAR 0.0693 0.0025 0.1312 0.4059 0.1386 0.0743 0.0569 0.0050 0.1163
PLUR 0.0534 0.0049 0.1019 0.2184 0.1748 0.0971 0.1359 0.0097 0.2039
IUT 0.1016 0.0070 0.1086 0.1384 0.1524 0.0946 0.2259 0.0140 0.1576

54
 Valeurs propres % d’inertie Cumul
1 0.040 0.837 0.837
2 0.005 0.115 0.952
3 0.001 0.024 0.976
4 0.001 0.020 0.996
5 0.000 0.003 0.999

attention à surinterpréter la distance entre les axes et se concentrer plutôt sur les angles
entre les points et les axes.

La dépendance entre le type d’études et l’origine sociale peut être résumée

par deux facteurs
• le premier oppose les études de médecine, caractéristiques des fils de
professions libérales et cadres supérieurs, aux études en IUT, caractéris-
tiques des fils d’ouvriers ;
• le deuxième oppose les fils d’exploitants agricoles à ceux de la CSP «
autres » et les études de pharmacie et d’IUT aux études de lettres.

Analyse des correspondances

0.4

EXPA
0.2

PHAR IUT
Facteur 2 (11.49%)

SCIE
MEDD OUVR
LICS PATR SEC
0.0

CADM SERV
DROI EMPL SALA
LETT PLUR
AUTR
-0.2

-0.4 -0.2 0.0 0.2 0.4 0.6

Facteur 1 (83.72%)

55
Tableau 8 – Cosinus carré des angles des vecteurs avec les sous-espaces

(a) Vecteurs lignes (b) Vecteurs colonnes

Facteur 1 Facteur 2 Facteur 1 Facteur 2

EXPA 0.50 0.47 DROI 0.00 0.31
SALA 0.93 0.01 SEC 0.03 0.03
PATR 0.06 0.15 LETT 0.59 0.40
LICS 0.99 0.01 SCIE 0.56 0.15
CADM 0.39 0.02 MEDD 0.98 0.00
EMPL 0.81 0.08 PHAR 0.61 0.32
OUVR 0.95 0.01 PLUR 0.85 0.11
SERV 0.81 0.00 IUT 0.87 0.11
AUTR 0.25 0.73

Tableau 9 – Contribution aux inerties associées aux axes factoriels

(a) Lignes (b) Colonnes

Facteur 1 Facteur 2 Facteur 1 Facteur 2

EXPA 0.07 0.48 DROI 0.00 0.03
SALA 0.01 0.00 SEC 0.00 0.00
PATR 0.00 0.02 LETT 0.06 0.31
LICS 0.53 0.03 SCIE 0.03 0.07
CADM 0.01 0.00 MEDD 0.50 0.02
EMPL 0.02 0.02 PHAR 0.06 0.24
OUVR 0.32 0.02 PLUR 0.02 0.01
SERV 0.02 0.00 IUT 0.33 0.31
AUTR 0.02 0.43

56
 Chapitre 6
Partitionnement de données
Dans les deux prochains chapitres, nous aborderons quelques techniques de partition-
nement de données et de discrimination. La partitionnement de données (“clustering”)
est une technique d’apprentissage non-supervisée qui consiste à grouper des observa-
tions. Les groupes sont constitués basés sur des mesures de dissimilarité (ou de dis-
tance) entre les individus. Les groupements sont inconnus, contrairement à l’analyse
discriminante que nous verrons dans le prochain chapitre.

À la base de la partitionnement de données, on trouve le concept de dissimilarité . Soit

xi , x j , xk ∈ R p . Une mesure de dissimilarité satisfait (1) d( xi , x j ) ≥ 0, (2) d( xi , xi ) = 0
et (3) d( xi , x j ) = d( x j , xi ). Si en plus d respecte l’inégalité du triangle, (4) d( xi , x j ) ≤
d( xi , xk ) + d( xk , x j ), alors d définit une distance.

On considère quelques distances, notamment

q
• Euclidienne : d( x, y) = ( x − y)> A( x − y) avec A = I ou plus généralement. On
prend souvent A = S−1
• Minkowski :
!1
p m
m
d( x, y) = ∑ | xi − yi | ,
i =1

appelée Manhattan (ou “city block”) quand m = 1 et euclidienne si m = 2.

p
• Canberra d( x, y) = ∑i=1 | xi − yi |/(| xi | + |yi |).

À partir de la distance, on peut construire facilement une similarité ,

1
sik = ,
1 + dik

où dik représente la distance entre l’objet i et l’objet k. Si on dispose d’une matrice de

similarité définie non-négative, on peut construire dik = 2(1 − sik ).
p

Un exemple de dissimilarité entre observations est d( xi , x j ) = 1 − ρij , où ρij est la cor-

rélation linéaire de Pearson.

6.1 Méthodes de regroupement hiérarchiques

Elles consistent à

• agglomérer des individus ou des groupes d’individus,

• partitionner de plus en plus.

Le résultat est présenté sous la forme d’un dendrogramme. Dans le dendrogramme

ci-dessus, on commence par obtenir la distance minimale entre les observations, ici

57
 F IGURE 5 – Exemple de dendrogramme

[1] [2] [3] [4] [5]

distance: échelle
de groupement

entre [1] et [2]. On continue itérativement, en groupant [3] et [4], puis le groupe [12] et
[34] et finalement en incluant [5] à [1234]. On dénote ci-dessous par dk(ij) la distance
entre l’observation k et les groupes i et j (formés respectivement de ni ≥ 1 et n j ≥ 1
observations).
Proposition 6.1 (Méthodes hiérarchiques)
1. distance unique (single linkage) : distance minimale dk(ij) = min(dki , dkj ), corres-
pond au plus proche voisin
2. distance complète (complete linkage) : distance maximale dk(ij) = max(dki , dkj ), voi-
sins les plus éloignés

F IGURE 6 – Illustration de la distance unique (bleu) et de la distance complète (rouge)

3. distance moyenne (average linkage) : moyenne de toutes les distances,

n1 n2
1
dk(ij) =
ni n j ∑ ∑ dij
i =1 j =1

4. barycentre (centroid) : distance entre les deux moyennes,

! 2
ni nj
dk(ij) = X̄ k − X̄ i + X̄
ni + n j ni + n j j

5. médiane
2
1
dk(ij) = X̄ k − ,


X̄ i + X̄ j
2

correspondant à la distance entre X̄ k et X̄ = ( X̄ i + X̄ j )/2.

58
 6. méthode de Ward : basée sur la somme des carrés (SC). La distance entre deux
groupements est la somme des carrés inter-groupes,

ni n j
dij = SCB = ni k X̄ i − X̄ k2 + n j k X̄ j − X̄ k2 = k X̄ i − X̄ j k2
ni + n j
avec
ni X̄ i + n j X̄ j
X̄ = .
ni + n j

Pour recalculer les distance entre les groupes, on utilise la formule de Lance et Williams,
soit

dk(i,j) = αi dki + α j dkj + βdij + γ|dki − dkj |.

Le tableau suivant résumé les relations pour les différents modèles.

αi αj β γ
1 1
Simple 2 2 0 − 12
1 1 1
Complet 2 2 0 2
ni nj
Moyenne ni + n j ni + n j 0 0
ni nj
Centroïde ni + n j ni + n j − αi α j 0
1 1
Médiane 2 2 − 14 0
n k + ni nk +n j
Ward n k + ni + n j n k + ni + n j − nk +nnki +n j 0

Dans ces méthodes hiérarchiques, le choix du nombre de groupes est souvent fait en
maximisant l’expression

tr( B
b ) / ( g − 1)
tr(Wc )/(n − g)

ou en examinant le dendrogramme et en choisissant un niveau adéquat. Cela permet

éventuellement de déterminer les aberrances si une observation rejoint le reste seule-
ment à la fin.

6.2 Méthodes non hiérarchiques

La problème général de la partitionnement de données optimale en un nombre prédé-

fini de g groupements est en général NP complet, puisque le nombre de partitions de
n individus en g groupes est de l’ordre O( gn /g!).

Une alternative consiste à procéder à une affectation initiale des observations en g grou-
pements et réassigner les individus de façon à optimiser un critère à chaque étape jus-
qu’à convergence du processus de réassignation. Le résultat dépend de la répartition

59
initiale des groupes, lequel n’est donc pas anodin, puisqu’il peut conduire à des solu-
tions finales différentes si l’algorithme se termine.
Les K-moyennes consiste à minimiser tr(W
c ) ou maximiser tr( B
b ).
Algorithme 6.1 (Partitionnement de données par K-moyennes (Lloyd))
Soit X un tableau d’observations n × p et K un nombre pré-spécifié de groupements.
1. Assignation initiale des observations aux K groupements et calcul des centroïdes
x̄k ou pré-spécification des centroïdes x̄k pour k = 1, . . . , K.
2. Calculer la distance Euclidienne carrée entre chaque observation et le centroïde x̄k
de son groupe d’appartenance c(xi ) ∈ {c1 , . . . , cK }

K
ESS = ∑ ∑ (xi − x̄k )> (xi − x̄k ).
k =1 i:c(xi )=ck

3. Assigner : ck ← {x j : kx j , x̄k k ≤ kx j , x̄i k, i = 1, . . . , K }, k = 1, . . . , K

4. Mise à jour des centroïdes : x¯k ← card−1 (ck ) ∑i:c(xi )=ck xi .
5. Itérer les étapes 3. et 4. tant que des observations sont réassignés à des groupes

L’algorithme de K-moyennes minimise le critère ESS, qui décroît avec chaque itéra-
tion et donc la convergence en temps fini est garantie. Comme indiqué précédemment,
l’obtention d’une solution optimale n’est pas garantie avec cet algorithme glouton et
dépend de l’allocation initiale des centroïdes ; en pratique, on essaiera plusieurs va-
leurs initiales. L’algorithme de K-moyenne tend à créer des groupes sphériques de taille
équivalente, puisque la structure de variance est la même pour chaque groupement.
On peut aussi procéder en utilisant d’autres critères de réallocation.
• minimiser |W
c |, maximiser | T
b | / |W
c|
−1
• maximiser tr(W
c B b)

Une variante des K-moyennes, les K-médianes, minimise plutôt la distance de Manhat-
tan (`1 ) entre les observations et la médiane de chaque groupe, à savoir |xk − medk (x)|.
Une autre variante utilise comme représentant de groupe une des observations, appelé
médoïde, en lieu et place des centroïdes. Le critère de minimisation est basé sur une
dissimilarité ou une distance. Ce critère est plus robuste face aux données aberrantes.
Deux variantes existent : le partitionnement autour de médoïdes (PAM) est le plus uti-
lisé, puisqu’il utilise un algorithme glouton.
Algorithme 6.2 (K-médoïdes et partitionnement autour de médoïdes)
Soit D = (dij ) une matrice n × n de proximité et K un nombre pré-spécifié de groupe-
ments.
1. Assignation initiale des items en K groupements et localisation du médoïde pour
chaque groupe. Le médoïde du groupe k est l’observation mk qui minimise la dissimi-
larité totale ∑i∈ck dimk .

60
 K-médoïdes :
(a) Pour le ke groupe, réassigner l’observation ik au groupe la plus près telle que la
fonction objective

K
ESSmed = ∑ ∑ diik
k =1 i:c(i )=ck

est réduite.
(b) Répéter l’étape 3 et la réassignation jusqu’à ce que plus aucune réassignation d’ob-
servations n’ait lieu.
Partitionnement autour de médoïdes :
(a) Pour chaque groupe, interchanger le médoïde avec l’observation qui réduit le plus
la fonction objective ESSmed
(b) Répéter les permutations jusqu’à ce que plus aucune réduction de ESSmed n’ait
lieu.

Un diagnostic graphique utile pour les méthodes de partitionnement de données est

la silhouette, qui est basée la matrice de proximité D. Pour chaque observation, la
moyenne des distances intra-classe est calculée pour l’observation i, dénotée ai , ainsi
que la plus petite distance moyenne inter-classe pour l’observation, dénotée bi . Plus
précisément, on a

1
ai =
card(c(i )) − 1 ∑ d(i, j)
j ∈ c (i )
1
bi = min
k∈{1...,K } card(ck ) ∑ d(i, j).
j∈ck
k 6 = c (i )

La largeur de la silhouette pour i pour une méthode de partitionnement de données

avec K groupements est donnée par

bi − a i
si,K =
max{ ai , bi }

et −1 ≤ si,K ≤ 1. Des grandes valeurs positives (ai ≈ 0) indiquent que la ie observation

est bien représentée par le groupe, des valeurs négatives (bi ≈ 0) indiquent une mau-
vaise allocation et si,K presque nulle indique que l’observation se situe entre deux grou-
pements. Une règle subjectives veut que si maxi si,K < 0.25 , aucune structure substan-
tielle n’a été détectée et qu’entre 0.26 et 0.50, la structure soit potentiellement artificielle.
Graphiquement, on recherche des groupes plutôt homogènes et une moyenne des sil-
houettes, n−1 ∑in=1 si,K , élevée.
Exemple 6.1 (Interprétation de silhouettes)
La figure 7 montre le résultat de la partitionnement de données à l’aide de l’algorithme

61
des K-moyennes sur un mélange de 4 lois binormales. Le diagnostic visuel indique
dans le cas K = 3 que certaines observations sont probablement plus susceptibles d’ap-
partenir à un autre groupe (en l’occurrence rouge). Le groupe vert est trop dispersé.
Avec K = 5, la séparation des groupements bleu-rouge est artificielle.

Longueur moyenne de silhouette: 0.72 K=3

5
0
x2

-5
-10

-10 -5 0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0

Longueur de silhouette si
x1

Longueur moyenne de silhouette: 0.86 K=4

5
0
x2

-5
-10

-10 -5 0 0.0 0.2 0.4 0.6 0.8 1.0

Longueur de silhouette si
x1

Longueur moyenne de silhouette: 0.77 K=5

5
0
x2

-5
-10

-10 -5 0 0.0 0.2 0.4 0.6 0.8 1.0

Longueur de silhouette si
x1

F IGURE 7 – Nuages de points et silhouettes pour la partitionnement de données de

données avec l’algorithme des K-moyennes

6.3 Méthodes basées sur un modèle

Soit un échantillon X = ( X 1> , . . . , X > >

n ) où X i est un p-vecteur colonne. On suppose
qu’on a K groupes, chacun caractérisé par une densité de dimension p, soit f k ( xi ; θk ) si

62
X i provient du groupe k pour k = 1, . . . , K. On réécrit la vraisemblance en fonction de
πk , la probabilité qu’une observation tombe dans le groupe k,

n K
L(θ1 , . . . , θK ; π1 , . . . , πK , X) = ∏ ∑ πk f k (xi , θk ).
i =1 k =1

Dans le cas de données manquantes (comme la classe d’une observation), on peut éli-
miner ces données, remplacer par la moyenne des colonnes (ce qui réduit la variance,
alors que les données manquantes devraient augmenter l’incertitude), faire une régres-
sion et prédire les données manquantes, voire faire de l’imputation multiple. La der-
nière solution consiste à utiliser l’algorithme de Dempster, Laird & Rubin (1977). Nous
utilisons l’algorithme d’espérance-maximisation puisque dans le cas présent l’informa-
tion est manquante.
On prend pour la densité

1 1
 
> −1
f k ( xi ; µ k , Σ k ) = exp − ( x − ) Σ ( x − ) .
2 i i
µ k k µ k
(2π ) p/2 |Σk |1/2

Nous allons prendre pour la matrice de covariance

Σk = λk Pk Ak P>
k (6.5)

avec Pk une matrice orthogonale de vecteurs propres, Ak une matrice diagonale dont
les éléments sont proportionnels aux valeurs propres et λk un scalaire. Pk gouverne
l’orientation de l’ellipsoïde, Ak sa forme et λk son volume, lequel est proportionnel à
p
λk det( Ak ).

Σk distribution volume forme orientation

λI sphérique constant constant non-définie
λk I sphérique variable constant non-définie
λPAP> elliptique constant constant constant
λk Pk Ak P> k elliptique variable variable variable
λPk AP> k elliptique constant constant variable
λk Pk AP> k elliptique variable constant variable

Le choix du nombre de composantes, de même que la forme du modèle, est basée sur la
sélection bayésienne de modèles et l’utilisation du facteur de Bayes. Spécifiquement, si
plusieurs modèles M1 , . . . , MO sont considérés avec chacun probabilité a priori p( Mo )
(souvent 1/O) pour un ensemble de données X, on a par le théorème de Bayes la loi
postérieure du modèle o ∈ {1, . . . , O}

p ( Mo | X ) ∝ p ( X | Mo ) p ( Mo )

63
où la vraisemblance intégrée p(X | Mo ) est obtenue en intégrant les paramètres incon-
nus du modèle et est donnée par la loi totale de probabilité
Z
p ( X | Mo ) = p(X | θo , Mo ) p(θo | Mo ) dθo

Le facteur de Bayes pour le modèle i versus j est Bij = p(X | Mi )/p(X | M j ) et la

comparaison est favorable au modèle i si Bij > 1 et il y a de très fortes preuves en
faveur de Mi si Bij > 100.

La vraisemblance intégrée p(X | Mo ) est approximée par le critère d’information bayé-

sien, ou BIC, donné par

BICo = −2 log p(X | Mo ) ≈ −2 log p(X | b

θo , Mo ) + νo log(n)

où νo = |θ| est le nombre de paramètres inconnus estimés dans le modèle o. On cherche

à minimiser le BIC. On peut aussi prendre comme critère de comparaisons des modèles
l’AIC ou le C p de Mallows.

Dans le cas de groupements ellipsoïdales, les mélanges de multinormales peuvent

mieux capturer certaines partitions inhomogènes, comme l’atteste cet exemple avec
ce mélange de normales en forme de souris. Le résultat de la partitionnement de don-
nées est indiqué par des cercles, et la couleur des points indique l’appartenance aux
différentes classes.

F IGURE 8 – Regroupemen pour le jeu de données mickey-mouse

Mélange de multinormales K-moyennes
0.8

0.8
0.6

0.6
0.4

0.4
0.2

0.2

0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8

64
6.4 Algorithme d’espérance-maximisation

L’algorithme d’espérance-maximisation est un algorithme itératif d’optimisation nu-

mérique, introduit par Dempster, Laird et Rubin (1977), qui permet d’obtenir les esti-
més du maximum de vraisemblance dans des problèmes avec données manquantes,
c’est-à-dire pour lesquelles les données du modèles sont partiellement observées. L’al-
gorithme produisant une séquence d’estimés qui converge vers l’estimé du maximum
de vraisemblance des données incomplètes à un rythme linéaire. Dans certains mo-
dèles, on peut aussi augmenter les données pour obtenir une vraisemblance complète
plus tractable et ainsi faciliter l’analyse.

Soit le couple (Y, U ) de données complètes. Les données U sont les données man-
quantes ou variables latentes, tandis que les données observées Y, aussi appelées don-
nées incomplètes. La vraisemblance pour les données complètes est

LY,U (θ; y, u) = f Y,U |θ (y, u; θ ),

tandis que la vraisemblance pour les données incomplètes est

Z
LY (θ; y) = f Y |θ (y; θ ) = f Y,U |θ (y, u; θ ) du.

L’algorithme EM facilite la maximisation de la vraisemblance pour les données incom-

plètes L(θ; Y ) en travaillant avec la vraisemblance des données complètes et la distri-
bution conditionnelle

f U,Y |θ (u, y | θ)
f U |Y,θ (u | y, θ) =
f Y |θ (y | θ)

On peut donc exprimer la log-vraisemblance des données incomplètes comme

`Y (θ; y) = `U,Y (θ; u, y) − `U |Y,θ (θ; u, y). (6.6)

Si on prend l’espérance de (6.6) par rapport à la densité conditionnelle f U |Y,θ (u | y, θ0 )

pour θ0 ⊂ Θ, on obtient
 
`Y (θ; y) = E fU |Y,θ `U,Y (θ; U, Y ) | y, θ0 − E fU |Y,θ `U |Y,θ0 (θ; U, Y ) | y, θ0


 
= Q(θ; θ0 ) − E fU |Y,θ `U |Y,θ0 (θ; U, Y ) | y, θ0

attendu que la log-vraisemblance des données observées ne dépend pas de U. On

conditionne sur une valeur spécifique θ0 lors du calcul de l’espérance conditionnelle.

L’algorithme d’espérance maximisation est un algorithme itératif qui produit une sé-
quence d’estimés qui convergent vers l’estimé du maximum de vraisemblance des don-

65
nées incomplètes. À partir d’une valeur initiale θ = θ0 , la r + 1e valeur de la séquence,
est construite selon
 
θ(r+1) = arg max E fU |Y,θ `U,Y (θ; U, Y ) | y, θ(r) = arg max Q(θ; θ(r) ).
θ∈Θ θ∈Θ

L’algorithme comporte deux étapes :

• étape E : calcul de l’espérance conditionelle de la log-vraisemblance pour les don-

nées complètes
• étape M : maximisation de ladite espérance.

Preuve On veut montrer que `Y (θ(r+1) ; y) ≥ `Y (θ(r) ; y). Puisque l’on choisit θ(r+1) qui
maximise la fonction Q, on a Q(θ(r+1) ; θ(r) ) ≥ Q(θ(r) ; θ(r) ) et il suffit de démontrer que
   
Ef `U |Y,θ(r) (θ(r) ; U, Y ) ≥ E f `U |Y,θ(r) (θ(r+1) ; U, Y )
U |Y,θ(r ) U |Y,θ(r )

Or, pour toutes fonctions de densité f 1 , f 2 pour une variable aléatoire U donnée, on a

f 2 (U )

 
E f1 (log f 1 (U )) − E f1 (log f 2 (U )) = −E f1 log
f 1 (U )
f 2 (U )
 
≥ − log E f1
f 1 (U )
f 2 (u)
Z
= − log f (u) du
f 1 (u) 1
Z
= − log f 2 (u) du = 0

par l’inégalité de Jensen. En appliquant ce résultat avec f 1 = f U |Y,θ(r) et f 2 = f U |Y,θ(r+1) ,

on obtient que la séquence {`Y (θ(r) )}r∞=1 converge vers un maximum local de la log-
vraisemblance `Y (θ) si la fonction est bornée. 

Des méthodes pour accélérer la convergence d’EM existent. Puisque l’algorithme EM

converge vers un maximum local, il peut être nécessaire de démarrer l’algorithme à
plusieurs valeurs de départ pour traiter les cas où la fonction objective est multimo-
dale. L’algorithme est très facile à coder, bien que les calculs pour l’espérance puissent
être difficiles en premier lieu. On peut déterminer la convergence en calculant (si dis-
ponible) |`Y (θ(r+1) ; y) − `Y (θ(r) ; y)| < ε ou |θ(r+1) − θ(r) | < δ pour des paramètres de
tolérance numérique δ et ε.
Exemple 6.2 (Données manquantes pour une variable binormale)
Soit W = (W1 , W2 )> ∼ N (µ, Σ) des données tirées d’une loi bivariée normale avec
moyenne µ = (µ1 , µ2 )> et matrice de covariance Σ = σσ11 σ12

21 σ22
. Les paramètres du
>
modèle sont θ = (µ1 , µ2 , σ11 , σ12 , σ22 ) . Supposons que l’échantillon obtenu contient
m données complètes w j = (w1j , w2j )> pour j = 1, . . . , m ainsi que w2j = (?, w2j )> ,
j = m + 1, . . . , m + m1 paire de données partiellement observées auxquelles manquent

66
la première coordonnée et finalement w1j = (w1j , ?)> , j = m + m1 + 1, . . . , n les m2 =
n − m − m1 paires auxquelles manquent la deuxième coordonnée.

La log-vraisemblance des données observées est donc

m
n+m 1 1
`(θ; W ) = −
2
log(2π ) − m log(det(Σ)) −
2 2 ∑ ( w j − µ ) > Σ −1 ( w j − µ )
j =1

1 1 n (w1j − µ1 )2 m+m1 (w2j − µ2 )2

− m1 log(σ22 ) − m2 log(σ11 ) − ∑ − ∑
2 2 j=m+m +1
2σ11 j = m +1
2σ22
1

Les données complètes consistent ici des n paires d’observations et la log-vraisemblance

complète est donnée par

n log(ξ ) 1
`c (θ; W, U ) = −n log(2π ) − − [σ22 T11 + σ11 T22
2 2ξ
− 2σ12 T12 − 2{ T1 (µ1 σ22 − µ2 σ12 ) + T2 (µ2 σ11 − µ1 σ12 )}
i
+n(µ21 σ22 + µ22 σ11 − 2µ1 µ2 σ12 )

où
n n
σ12
Ti = ∑ wij , Thi = ∑ whj wij , ξ = σ11 σ22 (1 − ρ2 ), ρ= √ .
j =1 j =1
σ11 σ22

où h, i = 1, 2. La distribution normale est une famille exponentielle régulière avec statis-

tique exhaustive ( T1 , T2 , T12 , T11 , T22 )> . Les estimateurs du maximum de vraisemblance
pour les données complètes sont

Ti T T T
bi =
µ , σhi = hi − h 2 i , h, i = 1, 2.
n n n
b

 
Pour calculer Q(θ; θ(r) ) = E log f (Y, U; θ) | y, θ(r) , on doit évaluer
 
2
Eθ(r) W1j | w2j , Eθ(r) W1j | w2j , j = m + 1, . . . , m + m1


 
2
Eθ(r) W2j | w1j , Eθ(r) W2j | w1j , j = m + m1 + 1, . . . , n

Or on sait par les propriétés de la loi multinormale que W2 | w1 suit une loi normale
 
−1
W2 | W1 = w1 ∼ N µ2 + σ12 σ11 (w1 − µ1 ), σ22 (1 − ρ2 ) .

Ainsi, l’étape E à l’itération r de l’algorithme EM prend la forme

(r )
(r ) (r ) σ (r )

Eθ(r) W2j | w1j = w2j = µ2 + 12 w1j − µ1


(r )
σ11

67
 (r )
2 (r )
 
2
Eθ(r) W2j | w1j = w2j + σ22.1 ,

pour j = m + m1 + 1, . . . , n. Le même raisonnement s’applique pour w1j pour j = m +

1, . . . , m + m1 .

Dans l’étape M de l’algorithme, on remplace les valeurs manquantes wij et wij2 pour
i = 1, 2 par leurs espérances conditionnelles, substituées dans les expressions pour
Ti , Thi , soit

(r ) (r ) (r ) (r )
(r +1) Ti (r +1) Thi T T
µi = , σhi = − h 2i , h, i = 1, 2.
n n n
On considère un cas simplifié où les observations manquent à un seul groupe.

Tableau 10 – Données manquantes pour une variable binormale)

w1 8 11 16 18 6 4 20 25 9 13
w2 10 14 16 15 20 4 18 22 ? ?

Tableau 11 – Premières 10 itérations pour l’algorithme EM de l’exemple 6.2

r µ1 µ2 σ11 σ22 σ12 −2`(θ(r) )

1 13 14.05537 40.2 130.53371 23.22854 125.27
2 13 14.47994 40.2 47.40122 21.20237 116.26
3 13 14.58502 40.2 30.89012 20.95190 113.49
4 13 14.60853 40.2 27.58288 20.89990 113.09
5 13 14.61375 40.2 26.92016 20.88843 113.03
6 13 14.61490 40.2 26.78735 20.88588 113.02
7 13 14.61516 40.2 26.76073 20.88532 113.01
8 13 14.61522 40.2 26.75539 20.88519 113.01
9 13 14.61523 40.2 26.75432 20.88516 113.01
10 13 14.61523 40.2 26.75411 20.88516 113.01

Exemple 6.3 (Mélanges finis)

Soit un échantillon X tiré d’un mélange fini de K composantes. La vraisemblance d’une
observation X i ∈ R p est alors

K K
f X i |θ ( xi | θ) = ∑ πk f k (xi | θk ), ∑ πk = 1
k =1 k =1

pour 0 < πk < 1. On peut augmenter la vraisemblance en considérant l’indicateur de

composante Ck , dénoté Zik = I (Ci = k ) pour chaque observation i et ainsi se débar-
rasser de la somme. Les données complètes sont donc constituées des tuples ( xi , zi )
avec zi = (zi1 , . . . , ziK ). Les zik correspond dans le cas présent aux variables latentes et

68
suivent une distribution multinomiale avec 1 essai. La log-vraisemblance peut s’écrire

n K
`(θk , πk ; z, x) = ∑ ∑ zik (log (πk ) + log ( f k (xi | θk ))) .
i =1 k =1

(r )
À l’étape r de l’algorithme, on impute zi par son espérance conditionnelle γi sachant
(r ) (r )
π (r) , θ1 , . . . , θK , qui découle de la formule de Bayes. La distribution conditionnelle
pour Ck est une distribution discrète sur {1, . . . , K } et puisque {Ci = k} ≡ { Zik = 1}

E f Z |X ,θ,π (I (Ci = k ) | xi , π, θ) = P f Z |X ,θ,π (Ci = k | xi , π, θ)

i i i i
que l’on peut écrire comme
πk f k ( xi | θk )
P f Z |X ,θ,π ( Zik = 1 | xi , π, θ) = K
= γk ( x i ) (6.7)
i i ∑ j =1 π j f j ( x i | θj )

où les variables latentes Zik sont conditionnellement indépendantes.

Dès lors,

Q(π, θ | π (r) , θ(r) ) = E f Z |X ,θ,π (`(θk , πk ; zik , x))

i i
K n K n
(r ) (r )
= ∑ ∑ γk (xi ) log(πk ) + ∑ ∑ γk ( xi ) log( f k ( xi | θk )) (6.8)
k =1 i =1 k =1 i =1

(r +1) (r +1)
On optimise alors chacune des termes séparément de façon à obtenir θ1 , . . . , θK
et π ( r + 1 ) . Le premier terme de l’équation (6.8) correspond à une vraisemblance multi-
nomiale, et donc comme précédemment, on obtient l’estimateur du maximum de vrai-
semblance
(r )
(r +1) ∑in=1 γk ( xi )
πk = (r )
.
∑Kj=1 ∑in=1 γ j ( xi )

Le deuxième terme peut être maximisé séparément, numériquement ou analytique-

ment selon la nature de f 1 , . . . , f K . On pose

n
(r +1) (r )
θk = arg max ∑ γk ( xi ) log( f k ( xi | θk ))
θk i =1

Dans le cas particulier où f k est une loi multinormale, les expressions pour le maximum
de vraisemblance de θk ≡ (µk , Σk ) pour k = 1, . . . , K sont explicites.
Algorithme 6.3 (Espérance-maximisation pour mélange fini de lois multinormales)
(0)
• Initialiser γk ( xi ) pour i = 1, . . . , n et k = 1, . . . , K
• Pour itération r = 1, . . . jusqu’à convergence
– Étape M : calculer les paramètres qui maximisent la vraisemblance étant donné

69
 (r −1)
γk ( xi ), soit

n (r ) n
(r ) (r −1) (r ) n (r ) 1 (r −1)
nk = ∑ γk ( x i ), πk = k ,
n
µk = (r ) ∑ γk ( xi ) xi
i =1 n k i =1

(r )
et Σk dépend de la structure de la matrice de covariance (voir équation (6.5)).
(r )
– Étape E : imputer les valeurs zik par γk ( xi ) étant donnés les paramètres estimés
à l’étape M,

(r ) (r ) (r )
π k f k µ k , Σ k | xi


(r )
γk ( x i ) = (r ) (r ) (r )
∑Kj=1 π j f j µ j , Σ j


| xi

Dans notre problème, il peut arriver que la matrice Σ b k soit mal conditionnée. 8 La
convergence de l’algorithme EM peut être lente si le nombre de groupe est élevé.

Estimation de la variance avec l’algorithme EM - cas régulier

Si les algorithmes de type quasi-Newton offrent un estimé de la matrice de variance

par le biais de la matrice hessienne, les écarts-types ne sont pas automatiquement dis-
ponibles avec l’algorithme EM. Nous assumerons subséquemment que des conditions
de régularité permettent d’interchanger dérivée et intégrale.

La fonction de score des données incomplètes est donnée par

 
∂ ∂
S(y; θ) = ` (θ; y) = E fU |Y,θ ` (θ; U, Y ) | y, θ
∂θ Y |θ ∂θ U,Y |θ
= E fU |Y,θ (Sc (u, y; θ)) (6.9)
 
∂ 0
= Q(θ | θ ) .
∂θ θ0 =θ

Cette autoconsistance de l’algorithme justifie le choix particulier de l’étape M ; au maxi-

mum de vraisemblance, l’estimé b θ résout l’équation de score S(y, θb) = 0.

On peut décomposer la variance en fonction de l’information pour les données com-

plètes et l’information manquante, données respectivement par

∂2
 
IU |Y (θ; y) = −E fU |Y,θ `U |Y (θ; U, Y ) | y
∂θ∂θ>
∂2
 
IU,Y (θ) = −E fU,Y |θ ` U,Y ( θ; U, Y ) .
∂θ∂θ>
bk = xl , Σ
8. S’il y a un groupe avec un seul individu xl , on aura alors µ b k → Od et la log-vraisemblance
tendra vers l’infini.

70
En dérivant équation (6.6) par rapport à θ deux fois, on obtient

∂2 ∂2 ∂2
` Y ( θ; y ) = ` U,Y ( θ; u, y ) − `U |Y,θ (θ; u, y).
∂θ∂θ> ∂θ∂θ> ∂θ∂θ>

En prenant l’espérance par rapport à f U |Y,θ à Y = y, on obtient

∂2
`Y (θ; y) = IU,Y (θ) − IU |Y (θ; y).
∂θ∂θ>

La première contribution est l’information maximale que l’on peut obtenir en n’ayant
pas observé U en prenant l’espérance sur tous les U possibles ; le deuxième terme est
la quantité d’information manquante, conséquence des données manquantes, que l’on
doit soustraire. En prenant l’espérance cette fois par rapport à f Y |θ, on obtient
 
IY (θ) = IU,Y (θ) − E fY |θ IU |Y,θ (θ; Y )

Un article de Louis (1982) montre que l’on peut exprimer l’information manquante
 
IU |Y,θ (θ; Y ) = E fU |Y,θ Sc (u, y; θ)Sc> (u, y; θ) − S(y; θ)S(y; θ)>

en utilisant équation (6.9). Or, la fonction de score au maximum de vraisemblance est

égale à 0, donc le terme de droite disparaît. On peut donc obtenir la matrice d’informa-
tion observée à bθ comme
 
I Y (b
θ) = IU,Y (b θ)Sc> (u, y; b
θ) − E f b Sc (u, y; b θ)
U |Y,θ
 
θ) .
θ) − Var f b Sc (u, y; b
= IU,Y (b U |Y,θ

Détection d’aberrances et diagnostics visuels

Une analyse exploratoire des données (à l’aide de nuages de points) peut être utile pour
diagnostiquer les groupements et détecter les aberrances. D’autres possibilités incluent
les profils (coordonnées parallèles), les diagrammes en étoiles, les faces de Chernoff,
les ACP et les diagrammes Andrews, correspondant à

x
f x (t) = √1 + x2 sin(t) + x3 cos(t) + · · ·
2

Les couleurs correspondent aux étiquettes des groupements. Si les courbes s’empilent,
alors les individus sont ressemblants.

71
F IGURE 9 – Diagrammes d’Andrews et en coordonnées parallèles pour le jeu de don-
nées iris
2
1
0
-1
-2

-3 -2 -1 0 1 2 3

Largeur sépales Largeur pétales Longueur sépales Longueur pétales

72
 Chapitre 7
Analyse discriminante et classification
L’analyse discriminante vise deux buts :
1. Séparation : décrire les caractères différentiels d’observations provenant de popu-
lations connues. Essayer de trouver les caractères discriminants qui séparent le plus
une collection d’individus.
2. Classification (ou allocation) : trier les observations en deux ou plusieurs groupes.
Obtenir une règle optimale d’allocation qui permet d’assigner un individu nouveau à
l’un des groupes.
En pratique, les deux buts se mélangent.

7.1 Le problème général de la classification

Soit deux populations Π1 , Π2 de densités f 1 ( x ), f 2 ( x ), respectivement et Ri la région
d’allocation à la population i. On dénotera par
Z
P ( R2 | Π1 ) = P ( X ∈ R2 | X ∼ Π1 ) = f 1 ( x ) dx
R2
Z
P ( R1 | Π2 ) = P ( X ∈ R1 | X ∼ Π2 ) = f 2 ( x ) dx
R1

où R2 = Ω \ R1 pour Ω := R1 t R2 . Cela revient à classifier un individu de la popu-

lation Π2 en 1 et réciproquement. On dénote également pi la probabilité a priori de
classer un individu dans Πi et C Ri | Π j le coût de mauvaise classification d’une ob-

servation de Π j dans Πi . Le coût espéré de mauvaise classification, dénoté CEMC est

F IGURE 10 – Deux régions d’allocation en 2D

R1

R2

donné par

CEMC = C ( R2 | Π1 ) P ( R2 | Π1 ) p1 + C ( R1 | Π2 ) P ( R1 | Π2 ) p2

73
Proposition 7.1
Les régions R1 et R2 qui minimisent ce coût CEMC sont

f1 (x) C ( R1 | Π2 ) p2 f1 (x) C ( R1 | Π2 ) p2
   
R1 := ≥ , R2 := < .
f2 (x) C ( R2 | Π1 ) p1 f2 (x) C ( R2 | Π1 ) p1

En pratique, on sélectionne souvent des probabilités à priori égales, de sorte qu’on a le

ratio des densités à comparer avec 1.

F IGURE 11 – Discrimination de deux populations et erreur de classification

R1 R2
0.8
Densité

Π1 Π2
0.4
0.0

-2 0 2 4 6

Preuve Par définition,

Z Z
CEMC = C ( R2 | Π1 ) p1 f 1 ( x ) dx + C ( R1 | Π2 ) p2 f 2 ( x ) dx
R2 R1
 Z  Z
= C ( R2 | Π1 ) p1 1 − f 1 ( x ) dx + C ( R1 | Π2 ) p2 f 2 ( x ) dx
R1 R1
Z
= (C ( R1 | Π2 ) p2 f 2 ( x ) − C ( R2 | Π1 ) p1 f 1 ( x )) dx + C ( R2 | Π1 ) p1
R1

en notant que tous les termes sont positifs. Le minimum sera donc atteint lorsque l’in-
tégrant est plus petite ou égale à zéro. Ainsi, R1 sera telle que C ( R1 | Π2 ) p2 f 2 ( x ) −
C ( R2 | Π1 ) p1 f 1 ( x ) ≤ 0, c’est-à-dire tel que

f1 (x) C ( R1 | Π2 ) p2
≥ .
f2 (x) C ( R2 | Π1 ) p1

Pour R2 , c’est évidemment le contraire. 

On peut adopter d’autres stratégies

74
 1. Minimiser la probabilité totale de mauvaise classification (revient au cas précédent
avec C ( R1 | Π2 ) = C ( R2 | Π1 ).
2. Allouer une nouvelle observation X0 à la population avec la plus grande probabilité
a posteriori

P ( X0 | Π 1 ) p 1 · p1 f 1 ( x0 )
P ( Π 1 | X0 ) = ∼
P ( X0 | Π 1 ) p 1 + P ( X0 | Π 2 ) p 2 p1 f 1 ( x0 ) + p2 f 2 ( x0 )
P ( X0 | Π 2 ) p 1 · p2 f 2 ( x0 )
P ( Π 2 | X0 ) = ∼
P ( X0 | Π 1 ) p 1 + P ( X0 | Π 2 ) p 2 p1 f 1 ( x0 ) + p2 f 2 ( x0 )

Cela revient au cas (1).

La généralisation à plus de deux populations (disons g) du CEMC allouera l’observa-

g
tion X0 = x0 à la population Πk pour laquelle ∑i=1,i6=k f i ( x ) pi C ( Rk | Πi ) est la plus
petite. Si C ( Rk | Πi ) = 1 (toutes sont égales), cette expression sera la plus petite si
pk f k ( x ) > pi f i ( x ) pour i 6= k.

Exemple 7.1 (Application du CEMC à deux lois normales)

Soit Π1 ∼ N (µ1 , Σ1 ) et Π2 ∼ N (µ2 , Σ2 )

1. Cas Σ1 = Σ2 = Σ : la densité est donnée par

1
 
f i ( x) = (2π )− p/2 |Σ|−1/2 exp − ( x − µi )> Σ−1 ( x − µi )
2
et donc
1 C ( R1 | Π2 ) p2
 
R1 : (µ1 − µ2 )> Σ−1 x − (µ1 − µ2 )> Σ−1 (µ1 + µ2 ) ≥ log
2 C ( R2 | Π1 ) p1
1 ( R1 | Π2 ) p2
 
C
R2 : (µ1 − µ2 )> Σ−1 x − (µ1 − µ2 )> Σ−1 (µ1 + µ2 ) < log
2 C ( R2 | Π1 ) p1

Dans ce premier cas, la discrimination est une droite formée par les points d’intersec-
tion des ellipses. En pratique, on utilisera les estimateurs usuels et on remplacera µ1
par x̄1 , µ2 par x̄2 et Σ par S p . La fonction de classification/discriminant d’Anderson
sera donnée par

1 C ( R1 | Π2 ) p2
 
1 1
( x̄1 − x̄2 )> S−
p x− ( x̄ − x̄2 )> S−
p ( x̄1 + x̄2 ) ≥ log
2 1 C ( R2 | Π1 ) p1

2. Le cas où Σ1 6= Σ2 mène à

1 C ( R1 | Π2 ) p2
 
R1 : − x> (Σ1−1 − Σ2−1 ) x + (µ1> Σ1−1 − µ2> Σ2−1 ) x − C ≥ log
2 C ( R2 | Π1 ) p1

75
 où C est donné par

1 | Σ | 1  > −1
log 1 +

C= µ1 Σ1 µ1 − µ2> Σ2−1 µ2
2 | Σ2 | 2

La structure de la courbe dépend de Σ1−1 − Σ2−1 , qui peut être une parabole ou une
ellipse, mais peut ne pas être définie positive.

Exemple 7.2 (Cas de plusieurs lois normales)

Supposons que C ( Rk | Πi ) = 1 ; allouer x à la population Πk si pk f k ( x ) > pi f i ( x) pour
i 6= k, c’est-à-dire allouer x à Πk pour laquelle

1 1
− log (|Si |) − ( x − x̄i )> Si−1 ( x − x̄i )> + log( pi ).
2 2
est le plus grand. Si on fait en plus l’hypothèse que les Σi sont égaux, x est alloué à la
population Πk pour laquelle

1 1
− ( x − x̄i )> S− >
p ( x − x̄i ) + log( pi )
2
est le plus grand. La quantité qui apparaît plus haut n’est rien d’autre que la distance
de Mahalanobis Di2 vue précédemment.

7.2 Analyse discriminante et classification par la méthode de Fisher

Nous entamerons l’étude de la méthode pour deux populations. Soit X = x un vecteur

d’observations de dimension p appartenant à une population, soit Π1 ou Π2 , tel que
E ( X | Π1 ) = µ1 et E ( X | Π2 ) = µ2 et Σ = Var ( X ) . L’idée est d’effectuer la transfor-
mation Y = a> X tel que µ1Y = a> µ1 , µ2Y = a> µ2 , σY2 = a> Σa qui séparera le mieux
les points projetés, c’est-à-dire celle qui maximisera

(µ1Y − µ2Y )2 a> (µ1 − µ2 )(µ1 − µ2 )> a ( a > δ )2

2
= = >
σY a Σa
> a Σa

où δ = µ1 − µ2 . On cherche donc à maximiser

( a > δ )2
max = δ> Σ−1 δ.
a a> Σa

En utilisant la décomposition de Choleski pour la matrice de covariance, Σ = LL> et

Σ−1 = L−> L−1 et en posant z = L> a, y = L−1 δ, on peut appliquer l’inégalité de
Cauchy–Schwartz (eq. 3.1) et ainsi obtenir

(a> δ)2 ≤ a> Σaδ> Σ−1 δ

76
avec égalité si a = kΣ−1 δ. Le maximum vaut donc δ> Σ−1 δ (d’où l’égalité) et il est
atteint en prenant a = Σ−1 δ. La combinaison linéaire y = (µ1 − µ2 )> Σ−1 x est appelée
fonction discriminante de Fisher.

La règle d’allocation est la suivante : posons

1 1
m= (µ1Y + µ2Y ) = (µ1 − µ2 )> Σ−1 (µ1 + µ2 )
2 2
et soit x à allouer à

Π1 si (µ1 − µ2 )> Σ−1 x ≥ m

Π2 si (µ1 − µ2 )> Σ−1 x < m

avec en pratique x̄1 , x̄2 et S p = ((n1 − 1)S1 + (n2 − 1)S2 )/(n1 + n2 − 2).

Les graphiques suivants illustrent la classification de 75 femmes en deux populations

(30 individus sains et 45 hémophiles) en fonction de l’activité et de la quantité de fac-
teur anti-hémophilique A (FAH) présent dans le sang, ainsi que la classification pour
un mélange de lois binormales.
0.6
0.4
x2 = log(antigènes de types FAH)

0.2
0.0
-0.2
-0.4

-0.6 -0.4 -0.2 0.0 0.2 0.4

x1 = log(activité FAH)

F IGURE 12 – Discrimination linéaire et méthode de Fisher pour données d’hémophilie

77
 4
2
0
y

-2
-4

-3 -2 -1 0 1 2 3 4

F IGURE 13 – Discrimination linéaire et méthode de Fisher pour trois populations binor-

males

En projetant dans une direction, on obtient le même critère dans le cas de lois multi-
normales en dimension deux. Le résultat diffère dans le cas général, comme l’illustre la
Figure 13.

Proposition 7.2 (Méthode de Fisher pour plusieurs populations)

L’idée est de déterminer un nombre restreint de combinaisons linéaires a1> X, a2> X, . . .
utiles pour séparer les populations Π1 , . . . , Π g pour g ≥ 2. On va supposer l’égalité des
matrices de variance-covariance Σ1 = · · · = Σ g = Σ. On dénotera

g
1
µ̄ :=
g ∑ µj
j =1
et
g
B= ∑ (µ j − µ̄)(µ j − µ̄)> ,
j =1

la matrice de variance inter-groupe. On considère Y = a> X, avec a> µ j = µ jY ainsi que

Var Y | Π j = a> Σa, µ̄Y = a> µ̄. Formons le quotient

>
g
2 g
 
∑ j=1 µ jY − µ̄Y a> ∑ j=1 µ j − µ̄ µ j − µ̄ a a> Ba
= = .
σY2 a> Σa a> Σa

78
Maximisons a> Ba sujet à la contrainte a> Σa = 1, avec Lagrange, soit
 
L = a> Ba − λ a> Σa − 1
et donc
∂L
= 2Ba − 2λΣa = 0,
∂a

qui donne Ba = λΣa ou Σ−1 Ba = λa. On retombe sur un problème de valeurs propres–
vecteurs propres. La maximisation successive (avec orthonormalité par rapport à Σ) de
cette expression est donnée par les vecteurs propres de Σ−1 B et est atteinte pour les
min( p, g − 1) vecteurspropres correspondants
 (car rang( B) = g − 1). On vérifie que
Var ai> X = 1 et Cov ai> X, a> .


j X = δij

En pratique, on calcule les valeurs propres et les vecteurs propres de Σ−1/2 BΣ−1/2 ,
puisque comme lors de l’analyse canonique, la matrice de départ n’est pas symétrique.
Si λi est valeur propre de Σ−1/2 BΣ−1/2 pour i = 1, . . . , k = min( g − 1, p) et ei est
vecteur propre de Σ−1/2 BΣ−1/2 , alors ai = Σ−1/2 ei et de ce fait ai> Σ1/2 Σ1/2 a j = δij .

Dans le cas d’un échantillon, on remplace µi par X̄ i et Σ par

∑ i ( n i − 1) S i Wc
Sp = =
∑i i ( n − 1 ) ∑i i − g
n
avec comme à l’accoutumée
g
∑ j=1 n j X̄ j
X̄ = g
∑ j =1 n j

>
∑ nj


b=
B X̄ j − X̄ X̄ j − X̄
j
g
c=
W ∑ ( n j − 1) S j .
j =1

La fonction Y i = b ai> x, où b
ai est le vecteur propre correspondant à la valeur propre b
λi
−1 b − 1 > e
de S p B ∝ W B avec b
c b a S pb
i a j = δij , est appelé i fonction discriminante.

Allocation

Considérons une observation hypothétique X et Yj = a> j X pour j = 1, . . . k = min( g −

1, p), la transformation correspondant à chacun des vecteurs propres a j . Alors Y =
(Y1 , . . . , Yk )> est, par construction, tel que Var (Y ) = Ik . Le vecteur Y a pour moyenne
µlY = (a1> µl , a2> µl , . . . , a> >
k µl ) si l’observation provient du groupe l. Il est alors appro-
prié de calculer l’éloignement de Y = y par rapport au centre de chacune des popula-
tions µlY , soit (y − µlY )> (y − µlY ). On attribuera Y = y à la population pour laquelle
cette distance est minimisée.

79
On peut aussi se résoudre à n’utiliser que r < k fonctions discriminantes, comme lors
d’une ACP.

7.3 Autres considérations sur la discrimination

1. L’analyse discriminante peut s’effectuer comme une analyse canonique particu-

lière. On écrit la matrice de données n × p comme X = (X1 , . . . , X g )> où Xi est mi × p
avec la matrice Y de dimension n × ( g − 1) donnée par

1m1 0m1 0m1

 
···
 0m2 1m2 · · · 0m2
 
Y=

 .. .. .. .. 
 . . . . 


0m g 0m g ··· 1m g

∗ ∗ −1 ∗
b et nS∗ S22
tel que yij = 1 si xi appartient au groupe j, et nS11 =T 12 S21 = B.b Aussi
−1 −1
T
b B b a valeur propre γbi , W
c B b a valeur propre b
λi et ces dernières sont reliées par la
relation λi = γ
b bi /(1 − γbi )
2. L’analyse discriminante est très fortement liée au test MANOVA (moyennes diffé-
rentes pour faire une discrimination).
3. Il existe une technique de régression pas à pas permettant de faire une analyse
discriminante.

7.4 Évaluation de la qualité de la discrimination

La validation croisée ou le calcul empirique du taux de classification sont des moyens

de déterminer la qualité de la règle.
1. Cas normal : On considère le cas de deux populations, avec Π1 : X ∼ N p (µ1 , Σ) et
Π2 : X ∼ N p (µ2 , Σ). On a

1
 
h( x) = a> x− ( µ + µ2 )
2 1
avec
1
 
P ( R1 | Π2 ) = P ( R2 | Π1 ) = Φ − ∆
2
où
∆ 2 = ( µ 1 − µ 2 ) > Σ −1 ( µ 1 − µ 2 )

avec les valeurs empiriques P b ( R1 | Π2 ) = P

b ( R2 | Π1 ) = Φ(−0.5∆ b 2 = ( X̄ 1 −
b ) avec ∆
1
X̄ 2 )> S−
p ( X̄ 1 − X̄ 2 ). En utilisant deux fois les données, on obtiendra un résultat (trop)
optimiste.
2. Resubstitution :
Soit n j la taille de Π j . On fait la discrimination, menant à l’allocation des données
initiales. Si nij dénote le nombre d’individus de Π j classés dans Ri (Πi ), avec comme

80
 estimé empirique P b = nij /n j . Même réserves pour cette méthode.
3. Validation croisée :
Effectuer la discrimination sur X −r , en ôtant l’observation numéro r, et calculer les
(r ) (r ) (r )
régions R1 , . . . , R g . Observer si l’observation r ∈ R1 . Soit ni1 le nombre d’indi-
(r )
vidus de Π1 pour lesquels xr ∈ Ri , i 6= 1, autrement dit sont mal classés. Alors
∗ = n∗ /n , pour i 6 = 1. Faire la même chose pour les échantillons 2, . . . , g. Cette
pi1 i1 1
méthode est meilleure au niveau de l’évaluation, mais est intensive en calculs.

7.5 Discrimination et régression logistique

Considérons la classification sur la base de variables explicatives X d’observations ap-

partenant à 2 groupes. On suppose que Y ∈ {0, 1} suit une loi Bernoulli, dont la densité
en caractérisation exponentielle est

p
   
f Y (y; p) = exp y log + log(1 − p)
1− p
 
= exp yθ + log(1 + eθ ) .

avec comme paramètre naturel θ = logit( p) := log( p) − log(1 − p). La moyenne de la

distribution est p = expit(θ ) := eθ /(1 + eθ ) ∈ (0, 1).

Un modèle généralisé linéaire est un modèle de régression de la forme

g(E (Y | X )) = η = α + X β, η∈R

pour une fonction lien g adéquate. On peut choisir par exemple la fonction lien na-
turelle g(·) = logit(·), qui correspond à p = expit(η ). Le modèle binomial résultant
satisfait les contraintes inhérente à l’hypothèse distributionnelle, dans le cas présent
E (Y | X ) ∈ (0, 1), tout en relâchant l’hypothèse de symétrie des résidus inhérente au
modèle linéaire.

La discrimination logistique est antérieure aux modèles généralisés linéaires et est

analogue à la discrimination linéaire. Soit la probabilité d’appartenance conditionnelle
à une population k ∈ {1, 2}, P ( Gk | X ). Par le théorème de Bayes, on peut écrire

P ( X | G1 ) P ( G1 )
P ( G1 | X ) =
P ( X | G1 ) P ( G1 ) + P ( X | G2 ) P ( G2 )
P ( X | G2 ) P ( G2 ) −1
 
= 1+
P ( X | G1 ) P ( G1 )

où P ( Gk ) dénote la probabilité a priori d’appartenir à la population k. On travaille alors

avec la vraisemblance des X. En supposant une discrimination logistique linéaire, en
assumant que les lois sont normales de même variance, soit X | G1 ∼ Nd (µ1 , Σ) et

81
X | G2 ∼ Nd (µ2 , Σ), on peut exprimer en forme close le maximum a posteriori, soit
 
P ( G1 | X ) = exp α + β> X P ( G2 | X )
  −1
P ( G2 | X ) = 1 + exp α + β> X

avec

β = Σ −1 ( µ 1 − µ 2 )
1 P ( G1 )
 
α = − (µ1 + µ2 )> Σ−1 (µ1 − µ2 ) + log
2 1 − P ( G1 )

Cette formulation nous force à imposer une hypothèse distributionnelle sur les X. La
régression logistique (avec η = logit( p) est plus parcimonieuse puisque l’on modélise
directement P Gj | X i et que seuls les coefficients de régressions sont spécifiés (la

variance d’une famille exponentielle étant fonction du paramètre naturel).

0.8
expit(β 0 + β 1 x )

0.4
0.0

-3 -2 -1 0 1 2 3

F IGURE 14 – Illustration d’un changement dans la probabilité lors d’un décalage du

taux de base (changement de α (rouge)) et de l’effet β (bleu) en fonction de la covariable
x.

On peut généraliser ce résultat au cas de K populations (données multinomiales) avec

une série de modèles

P ( Gk | X )
 
log = αk + β>
k X, k = 1, . . . , K − 1
P ( GK | X )

avec ici K comme classe de référence. Cela revient à imposer αK = 0, βK = 0 pour

satisfaire les contraintes d’identifiabilité ∑kK=1 P ( Gk | X ) = 1. La probabilité que l’ob-

82
servation i appartienne au groupement j est alors donnée par

exp(α j + β>
j Xi) exp(α j + β>
j Xi)
P Gj | X i = .


=
∑kK=1 exp(αk + βk X i )
>
1 + ∑kK=−11 exp(αk + β>
k Xi)

On attribue les nouvelles données X i au groupe maxk∈{1,...,K } P ( Gk | X i ).

7.6 Arbres de décisions

Un arbre de décision est une séquence de partitions binaires récursives qui mène au
partitionnement des données. On parle d’arbre de régression si la variable dépendante
est continue, et d’arbre de classification dans le cas de données catégoriques.

Arbres de classification

Soit Y ∈ {1, . . . , K } une variable catégorique et X un p-vecteur de variables expli-

catives. On considère une séquence de M règles de décisions binaires de la forme
X jm ≤ cm pour m = 1, . . . , M. Ces règles définissent des régions (des hyperrectangles)
dans l’espace des X dans lesquelles nous prédirons Y par un indicateur d’appartenance
à une catégorie k. Graphiquement, cette séquence peut être représentée par un arbre ;
chaque règle de décision crée un noeud τ scindant une branche en deux.

Arbre
Arbredede
Arbre de
classification
classification
classification desdes
des
iris iris
iris Classification des iris
7

versicolor
6
Longueur des pétales

long. pétales<2.45
≥2.45
4

larg. pétales<1.75
virginica
3

≥1.75

setosa versicolor virginica

2

50 0 0 0 49 5 0 1 45
setosa
1

0.5 1.0 1.5 2.0 2.5

Largeur des pétales

F IGURE 15 – Arbre de classification pour les données iris (gauche) et régions corres-
pondantes (droite).

Comment choisir les conditions booléennes pour séparer un noeud τm ? En sélection-

nant la variable X jm et une valeur cm de séparation, on partitionne la région Rm définie

83
à τm en deux hyperplans définis par

( g) (d)
Rm ( jm , cm ) = { X | X jm ≤ cm , X ∈ Rm }, Rm ( jm , cm ) = { X | X jm > cm , X ∈ Rm },

respectivement les branches droites et gauches du noeuds τm . La région Rm dépend

des ancêtres de τm .

Le choix optimal de la variable explicative et du point de séparation prend la forme

 

min min ∑ yi − I yi = k g + min ∑ |yi − I (yi = k d )| ,



jm ,cm kg ( g) kd (d)
xi ∈ Rm ( jm ,cm ) xi ∈ Rm ( jm ,cm )

un problème d’optimisation avec ici une fonction de coût binaire.

Si la région Rm compte Nm observations, la proportion d’observations du groupe k dans

ce noeud est donnée par

1
pbmk =
Nm ∑ I ( yi = k )
xi ∈ R m

et représente la probabilité conditionnelle qu’une observation tombe dans le groupe k

dans la région Rm . Si le coût de mauvaise classification est binaire, le choix optimal est
d’assigner les observations au groupe k représentant la classe majoritaire des observa-
tions dans la feuille m, soit au groupe k(m) = maxk pbmk = mink ∑ xi ∈ Rm |yi − 1yi =k |.
On peut ensuite énumérer toutes les séparations possibles (au plus nd choix) afin de
trouver les paramètres jm , cm optimaux et ainsi la règle de décision binaire optimale à
une étape donnée.

La qualité de la séparation cm au noeud τm dépend de la proportions d’observations

bien classées. Plusieurs mesures mesurant la qualité de la séparation existent : elles
sont appelées fonction d’impureté de noeud et dénotées Qm . Les trois principales sont

l’erreur de classification : 1 − pbmk ;

l’index de Gini : ∑kK=1 pbmk (1 − pbmk ) ;
l’entropie croisée : − ∑kK=1 pbmk log( pbmk ) (ou la déviance −2 ∑kK=1 pbmk log( pbmk )).

L’erreur de classification est moins sensible aux changements dans les probabilités d’as-
signation aux noeuds et est moins utilisée pour guider la croissance de l’arbre.

L’algorithme procède comme suit : à partir d’une première division, on traverse l’arbre
en séparant chaque noeud terminal intermédiaire jusqu’à saturer l’arbre, en divisant
une branche seulement si l’effectif de cette dernière est supérieure à un certain seuil
nmin .

On continue cette procédure jusqu’à saturer l’arbre, surajustant ce faisant le modèle. La

taille de l’arbre saturé T0 correspond au nombre de feuilles (noeuds terminaux) de T0 ,

84
dénotée | T0 |. Le choix optimal du sous-arbre de classification T ⊂ T0 est déterminé à
l’aide d’une critère prenant en compte la complexité du modèle et la qualité des pré-
dictions.

En indexant les noeuds terminaux de T par l = 1, . . . , | T | et les régions associées par

Rl , on définit le critère de coût complexité

|T |
Cα ( T ) := ∑ Nl Ql (T ) + α|T |,
l =1

où Ql ( T ) est une fonction d’impureté de noeud et α un paramètre de réglage qui dicte le

compromis entre taille de l’arbre et adéquation du modèle. Ce dernier est typiquement
choisi par validation croisée.

Pour chaque valeur de α, il existe un sous-arbre minimal unique, Tα , tel que Tα =

arg minT ⊂T0 Cα ( T ). Ce dernier peut être trouvé par une séquence d’itérations en fusion-
nant les noeuds intermédiaires et en choisissant celui qui produit la plus faible hausse
de ∑l Nl Ql ( T ), jusqu’à obtenir un arbre avec un seul noeud. Il est possible de démon-
trer que cette séquence de sous-arbres contient Tα . Les observations dans une feuille de
Tα seront attribuées au groupe présentant la plus forte proportion d’observations.

Forces et faiblesses des arbres de décision

1. Les solutions sont simples et faciles à implémenter. Il y a peu d’hypothèses proba-

bilistes. Par exemple la discrimination linéaire fait l’hypothèse d’égalité des matrices
de variance-covariance.
2. L’ajustement se fait par le biais d’un algorithme glouton (les décisions préalables ne
sont pas remises en cause et aucune recherche exhaustive de l’espace n’est effectuée).
3. La classification se fait variable par variable, ce qui est un avantage s’il y a des don-
nées manquantes. Si les variables sont fortement corrélées, les arbres de classifications
sont peu performants.
4. Les arbres de classification sont très instables et peu robustes et l’arbre n’a pas de
caractère inférentiel. On peut se protéger en le construisant avec un échantillon d’en-
traînement et en utilisant le reste des données pour estimer l’erreur.

Les arbres de décisions sont simplistes, mais sont à la base de plusieurs techniques mo-
dernes d’apprentissage statistique, comme les forêts aléatoires, les réseaux de neurones
et les machines à vecteurs de support.

85
 Chapitre 8
Copules
Dans ce dernier chapitre, nous aborderons la notion de copule, qui permet de regrouper
sous un même étendard les distributions multivariées.

8.1 Motivation et définition

Il existe peu de distributions multivariées canoniques, et la plupart de ces modèles

sont des généralisation de lois unidimensionnelles, dont le domaine et les marges sont
prescrits. Les exemples les plus connus sont les lois multinomiale (marge binomiale,
support Zd+ ), Dirichlet (marges Beta, avec pour support le simplexe unitaire Sd−1 ), Wi-
shart (marges χ2 , matrices positives semidéfinies) et elliptiques.

Deux questions motivent notre approche, soit

1. comment obtenir des fonctions de répartition multivariées valides avec des marges
arbitraires ?
2. comment comparer des données dont les marges diffèrent ?
Lemme 8.1
Soit F une fonction de répartition univariée et sa fonction de quantile,

F ← (u) := inf{ x : F ( x ) ≥ u}, u ∈ [0, 1].

Alors

1. F ← est non-décroissante et continue à gauche. Si F est continue, alors F ← est stric-

tement croissante (mais pas nécessairement continue).
2. F ← ◦ F ( x ) ≤ x, pour tout x ∈ R. Si F est strictement croissante, F ← ◦ F ( x ) = x.
3. F ◦ F ← (u) ≥ u pour tout u ∈ [0, 1]. Si F est continu, alors F ◦ F ← (u) = u

Preuve

2. F ← ◦ F ( x ) = inf{y ∈ R : F (y) ≥ F ( x )} =: A. Alors, x ∈ A et inf( A) ≤ x,

donc inf( A) = F ← ◦ F ( x ). Supposer F est strictement croissant, soit F ← ◦ F ( x ) < x
pour x ∈ R. Par contradiction ; en supposant qu’il existe a ∈ A tel que a < x. Mais
F ( a) < F ( x ) puisque F est strictement croissant, ce qui contredit le fait que a ∈ A et
donc F ( a) ≥ F ( x ).
3. F est par définition continue à droite. Nous avons égalité à 0, alors assumons sans
perte de généralité que y ∈ (0, 1] et F ← (1) < ∞. Alors F ← (u) ∈ R pour u ∈ (0, 1) et

F ( F ← (u)) = lim F ( x n ).
xn → F ← (u)
xn > F ← (u)

Maintenant xn > F ← (u), donc il existe xn∗ ∈ A := { x : F ( x ) ≥ u} tel que xn∗ < xn .

86
 F est non-décroissante, donc u ≤ F ( xn∗ ) ≤ F ( xn ) et de ce fait xn ∈ A, c’est-à-dire
F ( xn ) ≥ u. Ainsi F ◦ F ← (u) = limn→∞ F ( xn ) ≥ u. Si F est continue,

F ( F ← (u)) = lim F (yn ) ⇔ yn < inf( A)

n→∞
yn → F ← (u)
yn < F ← (u)

et donc yn ∈
/ A implique que F (yn ) < u.

Proposition 8.2 (transformée intégrale de probabilité)

Si Y est une variable aléatoire absolument continue de fonction de répartition G, alors

G (Y ) ∼ U (0, 1)

Proposition 8.3 (transformation des quantiles)

Si U ∼ U (0, 1) a une loi uniforme standard, alors

P ( G ← (U ) ≤ x ) = G ( x )

où l’inverse généralisée est définie comme F ← ( p) := min{ x : F ( x ) ≥ p}

Preuve de la proposition 8.2 Observer que pour tout u ∈ (0, 1), x ∈ R,

F ← (u) ≤ x ⇔ u ≤ F(x)

En effet, F est non-décroissante et donc si F ← (u) ≤ x, donc F ( F ← (u)) ≤ F ( x ) et par

notre lemme 8.1, u ≤ F ( F ← (u)), ce qui prouve la nécessité de la condition.

Pour la suffisance, si u ≤ F ( x ), alors x ∈ {y : F (y) ≥ u} et inf{y : F (y) ≥ u} ≤ x ⇒

F ← (u) ≤ x. Soit U ∼ U (0, 1) et fixons X = F ← (U );

F ( x ) = P ( X ≤ x ) = P ( F ← (U ) ≤ x ) = P (U ≤ F ( x )) = F ( x )

puisque F ( x ) a pour domaine [0, 1]. 

Preuve de la proposition 8.3 Nous avons

P ( F ← (U ) ≤ x ) = P (U ≤ F ( x )) = F ( x ), x ∈ R.

en utilisant le point 3 du lemme 8.1 

La transformation des quantiles sert à générer des variables pseudo-aléatoires.

87
Exemple 8.1 (Génération de variables exponentielles)
Si U ∼ U (0, 1), alors − log(U ) ∼ E (1).

Si les distributions marginales sont inconnues, elles devront être estimées. On peut
utiliser un estimateur nonparamétrique de F : la fonction de répartition empirique, Fen .
Ce choix est motivé par le résultat suivant :
Théorème 8.4 (Glivenko–Cantelli)
Soit X ∼ F et Fen ( x ) := n−1 ∑in=1 I ( Xi ≤ x )

a.s.
k Fen ( x ) − F ( x )k∞ −→ 0

Il n’est pas difficile de dériver la normalité asymptotique de Fen puisque les indicateurs
√ d
sont des variables Bernoulli et que donc n( Fen ( x ) − F ( x )) −
→ N (0, F ( x )(1 − F ( x )))
par application du théorème central limite.
Pour transformer non-paramétriquement les données à l’échelle (pseudo)-uniforme,
on calcule les rangs des observations, dénotés Ri que l’on renormalise par un facteur
(n + 1)−1 . Les données ainsi obtenues sont dénommées pseudo-observations. Le choix
de (n + 1)−1 plutôt que n−1 évite les problèmes à la bordure et n’a pas d’effet asymp-
totiquement. 9 On aura ainsi

e i = rang( Xi ) = Ri .
U
n+1 n+1

La transformée intégrale de probabilité sert dans les diagrammes P–P (voir défini-
tion 2.19) pour vérifier l’adéquation d’une hypothèse paramétrique. La fonction de
quantile est utilisée dans les diagrammes Q–Q.

Copule

Une copule lie une fonction de répartition et des lois marginales. Formellement, c’est
une fonction de répartition multivariée dont les marges sont uniformes.
Définition 8.5 (Copule)
Une copule est une application C : [0, 1]d → [0, 1] satisfaisant
1. C a des marges uniformes : C (u) = ui si u j = 1 ∀ j 6= i et ui ∈ [0, 1].
2. C (u) = 0 si ui = 0 pour au moins un i ;
3. C est ∆-monotone, c’est-à-dire est non-décroissante en toute composante,

1 1
∑dj=1 i j
P (U ∈ ( a, b]) = ∆(a,b] C := ∑ ··· ∑ (−1) C (v1i1 , . . . , vdid ) ≥ 0
i1 =0 i d =0

pour des vecteurs a, b ∈ [0, 1]d avec a ≤ b et v ji j = I i j = 0 a j + I i j = 1 b j .

9. En effet, si le support d’une variable aléatoire est R, alors Fe← (1) = ∞.

88
 La condition 3. équivaut à nécessiter que le volume entre deux points de Rd soit non-
négatif (voir définition 1.2). Si la densité de C (u) existe, cette condition revient à re-
quérir que c(u) ≥ 0 pour tout u ∈ (0, 1)d .

Exemple 8.2 (Copules fondamentales)

La copule antimonotone W(u1 , u2 ) = max{u1 + u2 − 1, 0} est une fonction de répar-
tition valide décrivant la dépendance négative parfaite (c’est-à-dire la fonction de ré-
partition de U, 1 − U). Inversement, la copule comonotone M := id=1 ui décrit la dé-
V

pendance parfaite, soit la fonction de répartition de (U, . . . , U ). Les deux copules sont
singulières. Un autre cas intéressant est celui de l’indépendance : les composantes d’un
vecteur aléatoire X sont mutuellement indépendantes (Xi ⊥ ⊥ X j pour tout i 6= j) si et
seulement si
d
Π (u) := C (u1 , . . . , ud ) = ∏ ui , u ∈ [0, 1]d
i =1

1 1 1

0.8 0.8 0.8

0.6 0.6 0.6

W(u1 , u2 )
M(u1 , u2 )

Π(u1 , u2 )

0.4 0.4 0.4

1 1 1
0.2 0.2 0.2
0.8 0.8 0.8
0 0 0
1 0.6 1 0.6 1 0.6
0.9 0.9 0.9
0.8 0.8 0.8
0.7 0.4 0.7 0.4 0.7 0.4
0.6 u1 0.6 u1 0.6 u1
0.5 0.5 0.5
0.4 0.2 0.4 0.2 0.4 0.2
u2 0.3 u2 0.3 u2 0.3
0.2 0.2 0.2
0.1 0.1 0.1
0 0 0 0 0 0

F IGURE 16 – Copules antimonotone (gauche), d’indépendance (centre) et comonotone

(droite)
1.0

1.0

1.0
0.8

0.8

0.8
0.6

0.6

0.6
U2

U2

U2
0.4

0.4

0.4
0.2

0.2

0.2
0.0

0.0

0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

U1 U1 U1

Exemple 8.3 (Copule de Marshall–Olkin)

La famille de copules bivariées de Marshall–Olkin est

 u 1− α1 u , si u1 1 ≥ u2α2 ,
α
1− α 2
C (u1 , u2 ; α1 , α2 ) = min(u1 1 u2 , u1 u12−α2 ) = 1
 u 1 u 1− α2 , si u1 1 ≤ u2α2 .
α
2

89
Ces copules sont à la fois absolument continues et dotées d’une singularité. Puisque

∂2  u − α1 , α
si u1 1 > u2α2
1
C ( u1 , u2 ; α1 , α2 ) =
∂u1 ∂u2  u − α2 , α
si u1 1 < u2α2 ;
2

la masse de la composante singulière est concentrée sur la courbe u1 1 = u2α2 ∈ [0, 1]2 .
α
1.0
0.8
0.6
u2

1
0.4

c(u1 , u2 )
0.2

0.5

1
0.0

0.8
0
0.0 0.2 0.4 0.6 0.8 1.0 1 0.6
0.8 0.4
0.6
u1 0.4 u1
0.2 0.2
u2 0 0

F IGURE 17 – Échantillon et graphique d’une copule de Marshall–Olkin avec

α1 = 0.8, α2 = 0.2, n = 1000

8.2 Théorème de Sklar

Pourquoi restreindre son attention aux distributions multivariées dont les marges sont
uniformes ? Le résultat suivant démontre qu’il existe une relation entre les fonctions de
répartitions F de marges F1 , . . . , Fd et les copules.

Théorème 8.6 (Sklar (1959))

Soit F une distribution conjointe avec des distributions marginales F1 , . . . , Fd . Alors, il
existe une copule C : [0, 1]d → [0, 1] telle que, pour tout x1 , . . . , xd dans R = [−∞, ∞],

F ( x1 , . . . , xd ) = C ( F1 ( x1 ), . . . , Fd ( xd )) (8.10)

Si les distributions marginales sont continues, alors C est unique. Autrement, C est
uniquement définie sur Im( F1 ) × · · · × Im( Fd ).

Inversement, si C est une copule et F1 , . . . , Fd sont des fonctions de répartition univa-

riées, alors la fonction définie par eq. (8.10) est une distribution conjointe de distribu-
tion marginales F1 , . . . , Fd .

90
Preuve Nous supposerons que les fonctions F1 , . . . , Fd sont continues. Le cas général
est traité dans Nelsen (2006).

Soit X ∼ F et Uj := Fj ( X j ) pour j ∈ {1, . . . , d}. Par la transformée intégrale de proba-

bilité, Uj ∼ U (0, 1) et de ce fait la fonction de répartition C est une copule. Puisque Fj
est croissante sur Im( X j ), X j = Fj← ( Fj ( X j )) = F ← (Uj ) pour j = 1, . . . , d. Ainsi, pour
x ∈ Rd .
 
F ( x) = P X j ≤ x j ∀ j = P Fj← (Uj ) ≤ x j ∀ j

= P Uj ≤ Fj ( x j ) ∀ j

= C ( F1 ( x1 ), . . . , Fd ( xd ))

C est donc une copule qui satisfait l’équation (8.10).

Soit U ∼ C et le vecteur X := ( F1← (U1 ), . . . , Fd← (Ud )). Alors, pour x ∈ Rd ,

 
P ( X ≤ x) = P Fj← (Uj ) ≤ x j ∀ j

= P Uj ≤ Fj ( x j ) ∀ j

= C ( F1 ( x1 ), . . . , Fd ( xd )).

F telle que définie par équation (8.10) est donc la fonction de répartition de X et ses
marges sont F1 , . . . , Fd par la transformation des quantiles. 

La continuité absolue est cruciale pour l’unicité de la copule.

Exemple 8.4 (Distribution bivariée avec marges Bernoulli)
Soit ( X1 , X2 ) suivant une loi bivariée de Bernoulli avec P ( X1 = k, X2 = l ) = 1/4 pour
k, l ∈ {0, 1}. On déduit que P X j = k = 1/2 pour k ∈ {0, 1} et l’image Im( Fj ) =

{0, 1/2, 1} pour j = 1, 2. Toute copule telle que C (1/2, 1/2) = 1/4 satisfait ce critère
(par exemple la copule d’indépendance C (u1 , u2 ) = u1 u2 ou bien encore la copule
diagonale C (u1 , u2 ) = min{u1 , u2 , (u21 + u22 )/2}.

Une propriété fondamentale des copules est leur invariance à des transformations stric-
tement croissantes des distributions marginales.
Proposition 8.7
Soit ( X1 , . . . , Xd ) un vecteur aléatoire dont les marges sont continues et une collec-
tion {φi }id=1 de fonctions strictement croissantes sur Im( Xi ), respectivement. Alors
(φ1 ( X1 ), . . . , φd ( Xd )) et ( X1 , . . . , Xd ) ont la même copule C.

Preuve Soit X ∼ F de distributions marginales F1 , . . . , Fd . Sans perte de généralité, il est

possible d’assumer que φj est continue à droite aux points dénombrables de continuité

91
(puisque X j a par hypothèse une fonction de répartition continue, φj ( X j ) change sur
un ensemble de mesure nulle). Puisque φj est croissant sur Im( X j ) et que X j a une
distribution continue, φj ( X j ) a également une distribution continue et

Fφj (X j ) ( x ) = P φj ( X j ) ≤ x


 
= P Tj ( X j ) < x = P X j < Tj← ( x )


 
= P X j ≤ φ← j (x)

= Fj (φ←
j ( x ))

pour x ∈ R. Cette dernière égalité implique que

 
P Fj (φ← X u j = P Fj ( X j ) ≤ u j ∀ j = C (u)


j ( φ j ( j ))) ≤ j ∀

où la dernière égalité découle de l’unicité de la copule associée à une fonction de répar-

tition F. 

Les implications du théorème de Sklar sont nombreuses. D’une part, il permet d’étu-
dier la dépendance en séparant les effets marginaux de la structure conjointe. En effet,

• les copules sont invariantes aux transformations marginales strictement monotones

croissantes sur l’image, puisqu’elles sont basées sur les rangs. Cela veut dire que l’on
peut bâtir un modèle avec des marges (continues) arbitraires.
• il pose le principe que toutes les distributions multivariées de dimension d dont
les marges sont absolument continues correspondent à une copule, la recherche de
modèles adéquats peut se réduire à celui d’une copule correspondante.
• les formules dérivées pour les copules sont valides pour tous les modèles peu im-
porte leurs marges, et le support borné facilite les dérivations analytiques.

D’autre part, le résultat bien que général n’est pas constructif. Mikosch (2006) articule
plusieurs critiques, dont nous retenons un échantillon sélectif :

• la classe des copules est trop large pour être utile en pratique.
• le choix d’une copule en est souvent un de convenance mathématique.
• les copules ne s’insèrent pas dans le cadre existant des processus stochastiques : ce
sont des modèles strictement statiques.

Proposition 8.8 (Propriétés des copules)

Si U ∼ C alors 1 − U ∼ C, la copule de survie de C. Pour d = 2,

C ( u1 , u2 ) = 1 − (1 − u1 ) − (1 − u d ) + C (1 − u1 , 1 − u2 )
= C (1 − u1 , 1 − u2 ) + u1 + u2 − 1.

92
 1.0
3.0
A
C

2.5

0.8
2.0

0.6
X2

X2
1.5

0.4
1.0

0.2
0.5

B
0.0

0.0
0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0

X1 X1

F IGURE 18 – Illustration de l’invariance des rangs à des changements de marges

Si C = C, C est une fonction symétrique radiale.

Le vecteur aléatoire X est dit échangeable si

d
( X1 , . . . , X d ) = ( X π ( 1 ) , . . . X π ( d ) )

pour toute permutation (π (1), . . . , π (d)) de (1, . . . , d).

Par extension, une copule C est échangeable si c’est la fonction de répartition d’un
vecteur U échangeable dont les marges sont uniformes, c’est à dire si C est symétrique.
Exemple 8.5 (Distribution de Marshall–Olkin)
Soit X1 , X2 deux variables aléatoires représentant la durée de vie de deux composantes.
Si des chocs surviennent suivant trois processus ponctuels de Poisson indépendants de
paramètres λ1 , λ2 , λ12 ≥ 0 (selon qu’ils affectent la composante 1, 2 ou les deux simul-
tanément). Le temps de ces chocs sont des variables exponentielles indépendantes de
paramètre λ1 , λ2 , λ12 respectivement et la fonction de survie de la paire ( X1 , X2 ) est
donnée par

F̄ ( x1 , x2 ) = P ( X1 > x1 , X2 > x2 ) = P ( Z1 > x1 ) P ( Z2 > x2 ) P ( Z12 > max{ x1 , x2 }) .

Les fonctions de survie univariées de X1 et X2 sont F̄1 ( x1 ) = exp(−(λ1 + λ12 ) x1 ) et

F̄2 ( x2 ) = exp(−(λ2 + λ12 ) x2 ).
En prenant α1 = λ12 /(λ1 + λ12 ) et α2 = λ12 /(λ2 + λ12 ), on trouve que la copule
1− α
de survie ( Z1 , Z2 ) est C (u1 , u2 ) = min(u1 1 u2 , u1 u12−α2 ). Ce modèle correspond à la
copule de Marshall–Olkin de l’exemple 8.3.

L’intervalle borné [0, 1]d est utile pour la visualisation et la distribution uniforme facilite
l’obtention des formules, puisque la forme de la densité est extrêmement simple. En

93
utilisant la transformation intégrale de probabilité et la fonction de quantiles, il est en
revanche possible d’obtenir des distributions avec des marges autres que uniformes.
Une telle construction est appelée méta copule..
3

3
2

2
1

1
0

0
-1

-1

-1
-2

-2

-2
-3

-3

-3
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
3

3
2

2
1

1
0

0
-1

-1

-1
-2

-2

-2
-3

-3

-3

-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
3

3
2

2
1

1
0

0
-1

-1

-1
-2

-2

-2
-3

-3

-3

-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3

F IGURE 19 – Méta copules (via Mvdc) avec lois marginales N (0, 1). De gauche à droite :
haut en bas : Plackett, Gaussienne, Student–t ; Gumbel–Hougaard, Clayton, Frank ; Ali-
Mikhail-Haq, Joe, Farlie-Gumbel-Morgenstern. Modèles (sauf FGM) calibrés pour que
τ ≈ 1/3.

En résumé : en variant Fi dans l’équation (8.10) on peut construire des distributions

dont les marges sont arbitraires. En variant C, on peut altérer la structure de dépen-
dance entre X1 , . . . , Xd .

94
8.3 Mesure de la dépendance

Le théorème de Sklar n’est pas un résultat constructif, et en pratique le choix d’une

copule est souvent question de jugement. Avant de présenter quelques modèles clas-
siques, nous révisons les mesures de dépendance qui puissent être utiles pour la cali-
bration d’un modèle aux estimés empiriques issus d’un échantillon.

Bornes de Fréchet–Hoeffding

Les distributions multivariées ont des bornes naturels en terme de dépendance. Dans
le cas bivarié, ces deux bornes correspondent à des copules et sont strictes. Dans le
cas d ≥ 3, la borne inférieure ne décrit plus un modèle valide. Les bornes de Fréchet–
Hoeffding, de la série 12, peuvent s’écrire pour toute copule C (u1 , . . . , ud )

d
( )
max ∑ ui + 1 − d, 0 ≤ C (u) ≤ min {u1 , . . . , ud } =: M(u)
i =1

et sont valides ponctuellement.

Mesures d’association

Comment mesurer le degré de dépendance entre deux variables X, Y, au vu des bornes

de Fréchet ? Clairement, la dépendance négative versus positive entraîne un choix de
modèle potentiellement différent. Avant d’aborder quelques mesures d’association, ré-
visons les propriétés désirables que ce dernier devrait satisfaire.
Propriétés désirables d’un coefficient d’association r selon Embrechts, McNeil, Strau-
mann (1999) :
P1. symétrique : r ( X, Y ) = r (Y, X )
P2. normalisée : −1 ≤ r ≤ 1
P3. r = 1 ⇔ comonotonicité, r = −1 ⇔ antimonotonicité
P4. Pour T : R → R strictement monotone sur Im( X ),

r ( X, Y ) si T est croissante
r ( T ( X ), Y ) =
−r ( X, Y ) si T est décroissante

P5. r = 0 si et seulement si X ⊥
⊥Y
Notez qu’il n’existe pas de mesure satisfaisant P4 et P5 simultanément. Cela est dû à
l’existence de distributions sphériques. On peut en revanche exiger que X ⊥ ⊥ Y im-
plique que r = 0. Ces conditions, ré-énoncées dans un ouvrage non-technique, suivent
les grandes lignes de Schweizer & Wolff (1981), qui proposent de baser une mesure
d’association sur la copule. Cela permet de respecter P4 (contrairement au ρ de Pear-
son) et l’estimé résultant est toujours bien défini puisque l’ensemble [0, 1] est borné (par

95
contraste, ρ nécessite des variances finies pour les lois marginales). Une ultime condi-
tion souhaitable, hormis celle qui ont déjà été énoncées, provient de Scarsini (1984) et
veut que
d
P6. si ( Xn , Yn ) −
→ ( X, Y ) alors r ( Xn , Yn ) → r ( X, Y )

Corrélation linéaire (ρ de Pearson)

Mesure de dépendance linéaire la plus connue, elle est adaptée aux lois elliptiques.
Facile à calculer, elle est directement interprétable pour la loi multinormale, pour la-
√
quelle ρ = σ12 / σ11 σ22 . Pour une distribution arbitraire, les propriétés suivantes sont
valides :

• les corrélations linéaires atteignables forment un ensemble fermé [ρmin , ρmax ] in-
cluant 0.
• la corrélation minimale (maximale) est atteinte dans le cas de variables antimono-
tones (comonotones).
d
• Cor ( X, Y ) = ±1 si et seulement si Y = ± aX + b pour a > 0, b ∈ R.
• X ⊥ ⊥ Y implique que Cor ( X, Y ) = 0. L’inverse est vrai si la loi est multinormale
(faux en général)

La corrélation linéaire n’est pas invariante aux transformations marginales strictement

monotones. Un dernier point notable est le fait que ρ n’est pas toujours bien défini :
l’existence de ρ( X, Y ) est acquise si Var ( X ) < ∞ et Var (Y ) < ∞.

Les quelques exemples suivants servent à démolir quelques mythes entourant la cor-
rélation linéaire.
Exemple 8.6 (Lois de Pareto)
Concernant les derniers points, si X, Y proviennent d’une loi Pareto de paramètre 3 et
sont indépendantes, alors ρ( X, Y ) = 0, mais ρ( X 2 , Y ) n’existe pas.
Exemple 8.7 (Lois lognormales)
Reprenons l’exemple des bornes de Fréchet de la série 12, cette fois en fonction de
lois marginales Xi ∼ LN (0, σi2 ), i ∈ {1, 2}. Nous avons déjà démontré que les bornes
±1 ne sont pas toujours atteignables, comme les graphiques suivants l’illustrent. Si
σ12 = 1, σ22 = 25 alors ρ ∈ [−0.000002824, 0.000419091] !

Exemple 8.8 (Corrélation nulle)

Cet exemple montre que la corrélation de 0 n’implique pas l’indépendance. Soit la co-
pule

1
  
C (u1 , u2 ) = u1 u2 1 − 2θ u1 − (u1 − 1)(u2 − 1)
2

On peut montrer aisément que les marges sont uniformes et que ρ = 0 pour tout
θ ∈ [−1, 1], bien que la copule ne soit pas l’indépendance. Plus généralement : toute

96
 0 1

ρmax
ρmin
−0.5 0.5

−1 0
4 4
4 4
2 2 2 2
σ2 σ1 σ2 σ1
0 0

F IGURE 20 – Bornes de Fréchet–Hoeffding pour le ρ de Pearson de l’exemple 8.7

distribution symétrique autour de 0 avec 2e moments finis fournit un exemple de ce

point. Si X ∼ N (0, 1), Y = X 2 , alors Cor ( X, Y ) = E X 3 − E ( X ) E X 2 = 0.

$ de Spearman

Cette mesure d’association est équivalente à la corrélation linéaire de Pearson, mais à

l’échelle uniforme. Pour deux variables aléatoires X1 , X2 , il est défini comme

$( X1 , X2 ) = ρ( F1 ( X1 ), F2 ( X2 )).

En fonction de la copule, le $ de Spearman est égal à

E (U1 U2 ) − E (U1 ) E (U2 )

$ ( X1 , X2 ) =
Var (U1 ) Var (U2 )
p
Z 1Z 1
= 12 u1 u2 dC (u1 , u2 ) − 3
0 0
Z 1Z 1
= 12 C (u1 , u2 ) du1 du2 − 3
0 0

qui découle de l’identité de Hoeffding, c’est-à-dire

Z ∞ Z ∞
Cov ( X1 , X2 ) = ( F ( x1 , x2 ) − F1 ( x1 ) F2 ( x2 )) dx1 dx2
−∞ −∞

et du fait que E (U ) = 1/2, Var (U ) = 1/12 pour U ∼ U (0, 1).

Un estimateur empirique du coefficient de rang de Spearman est donné comme suit :

pour {rang( Xt,i ), rang( Xt,j )}, les colonnes i, j de X et t = 1, . . . , n,

n
12 1 1
    
$ ( Xi , X j ) =
n ( n2 − 1) ∑ rang Xt,i − (n + 1) rang Xt,j − (n + 1)
2 2
.
t =1

97
τ de Kendall

C’est une mesure nonparamétrique basée sur les rangs, qui mesure la concordance. Les
paires ( Xt1 , Xt2 ) et ( Xs1 , Xs2 ) sont concordantes si ( Xt1 − Xs1 )( Xt2 − Xs2 ) > 0.

Paire concordante Paire discordante

1.0

1.0
( Xt1 , Xt2 )
0.8

0.8
( Xt1 , Xt2 )
0.6

0.6
( Xs1 , Xs2 )
0.4

0.4
( Xs1 , Xs2 )
0.2

0.2
0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0

Définition 8.9 (τ de Kendall)

Le τ de Kendall est définie comme étant τ := P (concordance) − P (discordance), ou
ZZ
τ := 4 C (u1 , u2 ) dC (u1 , u2 ) − 1

Un estimateur empirique du coefficient de Kendall pour les colonnes i, j de X est

  −1
n
τn (Xi , X j ) = ∑ sign ( Xt,i − Xs,i )( Xt,j − Xs,j )


2 1≤ t < s ≤ n
Pn − Qn 4Pn
= n = −1
(2) n ( n − 1)

où Pn dénote le nombre de paires concordantes et Qn = n − Pn le nombre de paires

discordantes. La collection résulte en des matrices positives semi-définies.

Le τ de Kendall et le $ de Spearman sont basés sur les rangs et satisfont P1, P2 et P4.
En revanche, leur dérivation peut être ardue.

Note
Il est impossible de résumer la dépendance et de condenser toute l’information dans
un scalaire : les méta copules de la Figure 19 ont toutes le même τ de Kendall (à l’ex-
ception de la copule de FGM dont la borne supérieure est plus faible que τ = 1/3). La
plupart des mesures ne capturent pas les relations nonlinéaires complexes, comme un
nuage de point en étoile tel qu’illustré dans la figure 21. C’est pourquoi la représenta-
tion graphique des données demeure cruciale.

98
Coefficient de dépendance de queue

Il est possible de calculer la probabilité conditionnelle d’être dans le coin inférieur

gauche (supérieur droit) du carré unitaire, cette probabilité mesure la dépendance entre
les valeurs extrêmes. Dans le cas bivarié, le coefficient de dépendance de queue infé-
rieure est défini comme la limite

λl = lim P ( X2 ≤ F2← (q) | X1 ≤ F1← (q))

q →0+

P X2 ≤ F2← (q), X1 ≤ F1← (q)

= lim
P X1 ≤ F1← (q)


q →0+
C (q, q)
= lim
q →0+ q

si la limite λl ∈ [0, 1] existe. Similairement, le coefficient de dépendance de queue su-

périeure est

λu = lim P ( X2 > F2← (q) | X1 > F1← (q))

q →1−
1 − 2q + C (q, q)
= lim .
q →1− 1−q

Ce dernier est non-nul pour les modèles de valeurs extrêmes.

8.4 Familles et modèles

Les classes de modèles de copules les plus communes sont

1. induites ou implicites (copule découlant du théorème de Sklar, comme par exemple

les copules elliptiques ou de valeurs extrêmes)
2. explicites (construction directe de C ; notamment les copules archimédiennes).

Nous revoyons à tour de rôle les différentes classes de modèles et les propriétés défi-
nissant ces dernières.

8.4.1. Copules elliptiques

Bien connues, ces copules sont dites induites, c’est-à-dire que leur construction découle
du théorème de Sklar et des modèles définis précédemment (définition 2.17). Quelques
métacopules gaussiennes (modèles dont les marginales sont transformées d’uniformes
à normales) sont présentées en figure 21.

Les copules elliptiques sont populaires pour plusieurs raisons. D’abord, leur densité est
explicite, bien que la copule soit implicite (sous forme d’intégrale, comme la fonction
de répartition de la loi normale). Par exemple, la fonction de répartition de la copule

99
 Normale Student-t, ν = 1 Student-t, ν = 3
3

3
2

2
1

1
0

0
-1

-1

-1
-2

-2

-2
-3

-3

-3
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3

F IGURE 21 – Contour de la densité des métacopules pour des lois elliptiques (marges
N (0, 1)).

gaussienne pour d = 2
Z Φ −1 ( u ) Z Φ −1 ( u 2 )
1 x2 + y2 − 2ρxy
 
1
CN (u; ρ) = p exp − dy dx
1 − ρ2 −∞ −∞ 2(1 − ρ2 )

tandis que la fonction de densité est

1 2ρΦ−1 (u1 )Φ−1 (u2 ) − ρ2 (Φ−1 (u1 )2 + Φ−1 (u2 )2 )

 
cN (u; ρ) = p exp .
1 − ρ2 2(1 − ρ2 )

Les copules elliptiques sont faciles à échantillonner et sont utilisées couramment en

pratique de ce fait. Une restriction est la symétrie radiale qui associe la fonction de
répartition centrée à sa fonction de survie, C̄ (u) = C (u) ; cette propriété a des impli-
cations sur les coefficients de queue, qui doivent être égaux (λl = λu ). Des exemples
connus sont les copules gaussienne, tν , normale inverse gaussienne (NIG), hyperbo-
lique symétrique généralisée, Laplace, etc.

Exemple 8.9 (Copules elliptiques)

On peut considérer des distributions où µ = 0d et Σ = ℘(Σ) où ℘(Σ) est la matrice de
corrélation. La densité est alors

1 1 >
 
−1
f ℘(Σ) = g − x (℘(Σ)) x .
|℘(Σ)|1/2 2

Les dérivations peuvent découler de la forme quadratique R2 = X > (℘(Σ))−1 X ; voir le

tableau 1. Le coefficient de dépendance de queue de la copule gaussienne est 0 (ce qui
en fait un très mauvais choix pour la modélisation de retours financiers), tandis que

100
pour la copule de Student-t,
s
√ 1 − ρij
λl = λu = 2tν+1 ν + 1
1 + ρij

où ρij est l’élément hors-diagonale de ℘(Σ) (voir Demarta & McNeil (2005)).

Normal, ρ = 0.5 Student-t, ν = 2, ρ = 0.5 Cauchy, ρ = 0

1.0
1.0

1.0
0.8
0.8

0.8
0.6
0.6

0.6
U2

U2

U2
0.4

0.4
0.4

0.2

0.2
0.2

0.0
0.0
0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

U1 U1 U1

F IGURE 22 – Échantillons de copules elliptiques

Proposition 8.10 (Relation entre mesures de concordance)

Nous notons au passage quelques identités fort utiles : une relation entre ρ et τ, valide
pour toute copule elliptique, de même qu’une relation entre ρ et $ qui elle n’est valide
que pour la copule gaussienne. Elles sont données respectivement par

2
τ= arcsin(ρ)
π
6 ρ
$ = arcsin
π 2

8.4.2. Copules archimédiennes

Définition 8.11 (copule archimédienne)

Une copule copule archimédienne admet la représentation

Cψ (u) = ψ(ψ← (u1 ) + · · · + ψ← (ud )) (8.11)

où ψ : [0, ∞) → [0, 1], le générateur archimédien de la copule, est une application d-

monotone avec ψ(0) = 1 et limx→∞ ψ( x ) = 0 (si strict). Voir McNeil & Nešlehová (2009)
pour plus à ce sujet et pour une représentation stochastique.

101
 1.0
ρ
τ
$

0.5
0.0
r

-0.5
-1.0

-1.0 -0.5 0.0 0.5 1.0

F IGURE 23 – Relation entre ρ, $ et τ pour la copule gaussienne

Le (pseudo)-inverse ψ← : [0, 1] → [0, ∞] est donné par l’inverse de ψ sur (0, 1] et par

ψ← (0) = inf{ x : ψ( x ) = 0}.

Dans le cas bivarié, on a

Définition 8.12 (copule archimédienne bivariée)

Cψ (u1 , u2 ) = ψ(ψ← (u1 ) + ψ← (u2 ))

pour une fonction convexe ψ telle que ψ(0) = 1 et limx→∞ ψ( x ) = 0.

Les copules archimédiennes sont très populaires, notamment parce qu’elles sont ap-
parues très tôt dans la littérature. Elles ont souvent un ou deux paramètres de dépen-
dance, ce qui facilite l’estimation mais limite leur applicabilité. De nombreuses pro-
priétés extrémales et de dépendance sont disponibles à partir du générateur ψ. Elles
sont échangeables (symétriques). Certaines familles couramment utilisées incluent les
copules de Gumbel–Hougaard, Joe, Clayton, Ali–Mikhail–Haq et Frank. Pour les cinq
familles, les densités sont disponibles sous forme fermée et l’erreur moyenne quadra-
tique est ∝ 1/nd si les marges sont connues, voir l’article de Höfert, Machler et McNeil
(2012) pour des indications à ce sujet et les difficultés numériques sous-jacentes. Le
modèle est en revanche difficilement réaliste en dimension élevée, mais des modèles
de copules imbriquées ont commencé à émerger dans la littérature.

Une formule pour le τ de Kendall est disponible, dérivée par Genest & Rivest (1993),

102
est
Z 1
ψ← (t)
τ = 1+4 dt
0 (ψ← (t))0

tandis que les coefficient de dépendance de queue sont donnés par

ψ0 (2q)
λu = 2 − 2 lim
q →0+ ψ0 (q)
ψ0 (2q)
λl = 2 lim
q→∞ ψ0 (q)

Tableau 12 – Générateurs et propriétés de quelques copules archimédiennes

Famille ΘC ψ(t)
A [0, 1) (1 − θ ) / ( e t − θ )
C (0, ∞) (1 + θt)−1/θ
F (0, ∞) −log(1 − (1 − −θ −t
− 1/θ

e )e )/θ
G [1, ∞) exp − t
J [1, ∞) 1 − (1 − e−t )1/θ

Famille τ Im(τ) λL λU
1
A 1 − [2( θ + (1 − θ )2 log(1 − θ )]/3θ 2 (0, 3) 0 0
C θ/(θ + 2) (0,1) 2−1/θ 0
Rθ
F 1 − 4 0 t/[θ (et − 1)] dt − 1 /θ (0,1) 0 0

G (θ − 1)/θ [0,1] 0 2 − 21/θ

J 1 − 4 ∑∞k=1 1/ [ k ( θk + 2)( θ ( k − 1) + 2)] [0,1] 0 2 − 21/θ
Les familles sont celles de Ali-Mikhail-Haq (A), Clayton (C), Frank (F), Gumbel-Hougaard (G) et Joe (J).

Nous nous cantonnerons dorénavant à l’étude de fonctions dites complètement mono-

tone et sont valide pour tout d ∈ 2, 3, . . ., comme celles présentées dans le tableau 12.
Définition 8.13 (fonction complètement monotone)
Une fonction est complètement monotone si (−1)k ψ(k) ( x ) ≥ 0 pour tout k ∈ Z+ et ψ
est continue sur [0, ∞].

Par le théorème de Bernstein, la transformée de Laplace d’une variable aléatoire posi-

tive est complètement monotone, et donc ψ = L( F ) est un candidat pour un généra-
teur archimédien. Pour l’échantillonnage, on procède généralement en utilisant l’algo-
rithme de Marshall & Olkin (1988).
Algorithme 8.1 (Algorithme de Marshall-Olkin (1988))
Si ψ est complètement monotone :

1. Échantillonner X0 ∼ F0 = L−1 (ψ)

103
 2. Échantillonner X, où X j ∼ U (0, 1) pour j = 1, . . . , d.
3. Retourner Ui = ψ(− log( X j )/X0 ) pour j = 1, . . . , d.

8.4.3. Copules de valeurs extrêmes

En dimension d = 2, ces copules sont de la forme

log(u2 )
  
C A (u1 , u2 ) = exp log(u1 u2 )A
log(u1 u2 )

où la fonction de Pickands A : [0, 1] → [1/2, 1] est une application convexe telle que

max(t, 1 − t) ≤ A(t) ≤ 1, t ∈ [0, 1]

Comme leur nom l’indique, les copules de valeurs extrêmes ont des connections avec
les distributions max-stables, auxquelles elles sont complètement équivalentes.

8.4.4. Exemples de modèles archimédiens

Définition 8.14 (Copule de Gumbel–Hougaard)
C’est une copule archimédienne et de valeur extrême (correspondant à la distribution
de valeur extrêmes logistique), dont le générateur est donné par

ψ(t) = exp(−t1/θ ), ψ← (t) = (− log(t))θ

pour θ ∈ [1, ∞) et la copule s’écrit comme

 !1/θ 
d
Cθ (u) = exp − ∑ (− log(ui ))θ .
i =d

Cas limites :

• CGu (u | θ = 1) = Π (u), la copule d’indépendance et

• limθ →∞ CGu (u | θ ) = M(u).

τGu = (θ − 1)/θ et seule la dépendance positive est modélisée : τ ∈ (0, 1)

Définition 8.15 (Copule de Clayton)
Pour θ ≥ −1/(d − 1), le générateur

1/θ t−θ − 1
ψ(t) = (1 + θt)−
+ , ψ← (t) =
θ

strict si θ > 0, peut être vu comme transformée de Laplace d’une loi Γ(1/θ, 1/θ ). Le τ
de Kendall est τCl = θ/(θ + 2) et est borné entre −1/(2d − 3) et 1. D’autre propriétés
sont dérivées dans la série 13.

104
Note
La paramétrisation de la copule de Clayton présenté ci-dessus diffère de celle de R, qui
utilise plutôt ψ(t) = (1 + t)−1/θ comme générateur.

Clayton, τ = −0.5 Clayton, τ = 0.5

1.0

1.0
0.8

0.8
0.6

0.6
U2

U2
0.4

0.4
0.2

0.2
0.0

0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

U1 U1

Frank, τ = 0.5 Gumbel–Hougaard, τ = 0.5

1.0

1.0
0.8

0.8
0.6

0.6
U2

U2
0.4

0.4
0.2
0.2

0.0
0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

U1 U1

F IGURE 24 – Échantillons de copules archimédiennes

8.5 Estimation et inférence

Inférence nonparamétrique

Même l’estimation nonparamétrique copule en utilisant une extension multivariée de

la fonction de répartition empirique, il est préférable de spécifier un modèle paramé-
trique pour la copule. L’approche non-paramétrique repose sur le processus empirique
√
n(Cn − C ) ; ce dernier converge en loi vers un champ gaussien sous quelques hypo-

105
thèses de régularité, mais nécessite n grand. La copule empirique est néanmoins au
coeur de la plupart des tests d’indépendance et d’adéquation.

Définition 8.16 (copule empirique)

La copule empirique, proposée par Deheuvels (1979), est
 
n d Ri,j
1
Cn (u) = ∑I ≤ uj 
[
n n+1

i =1 j =1

Cette copule est l’équivalent de la fonction empirique en dimension d avec des marges
uniformes.

Inférence paramétrique

Soit X un échantillon n × d de vecteurs aléatoires X de fonctions de répartition F, ou de

façon équivalente de lois marginales continues F1 , . . . , Fd et de copule C. On suppose
que les données sont générées d’un modèle paramétrique F (θ0 ), où le vecteur de para-
mètres peut être partitionné en θ0 = (θ0,1 , . . . , θ0,d , θ0,C ) ∈ θ, avec Fj = Fj (·; θ0,j ) pour
j = 1, . . . , d et C = C (·; θ0,C ) et tel que Fj (; θj ) est continue pour tout θj ∈ Θ j .

Définition 8.17 (Densité d’une copule)

En vertu du théorème de Sklar, la densité de F (si elle existe) prend la forme

d
f ( x; θ0 ) = c( F1 ( x1 ; θ0,1 , . . . Fd ( xd ; θ0,d ); θ0,C ) ∏ f j ( x j ; θ0,j ),
j =1
où
∂d
c(u) = C ( u1 , . . . , u d ), u ∈ [0, 1]d
∂u1 · · · ∂ud

La log-vraisemblance basée sur X est

n
`(θ; X) = ∑ `(θ; X i )
i =1
n n d
= ∑ `C (θC ; F1 (Xi1 ; θ1 ), . . . , Fd (Xid ; θd )) + ∑ ∑ ` j (θj ; Xij ),
i =1 i =1 j =1

où

`C (θC ; u1 , . . . , ud ) = log c(u1 , . . . , ud ; θC )

` j (θj ; x ) = log f j ( x; θj ), j = 1, . . . , d.

106
L’estimateur du maximum de vraisemblance est donné par

θ = arg max `(θ; X)

b
θ∈Θ

Les modèles multivariés peuvent impliquer un grand nombre de paramètres si on op-

timise simultanément les distributions conjointes et marginales. L’optimisation multi-
variée peut être également numériquement instable.

Pour cette raison, la plupart des praticiens procèdent en deux étapes : estimation des
lois marginales d’abord, puis des paramètres de la copule.

Joe et Xu (1996) ont proposé une approche pour l’inférence en deux étapes.

Proposition 8.18 (Inférence pour les marges)

1. Estimer θ0,j par son maximum de vraisemblance b
θj à partir de la log-vraisemblance
n
marginale ∑i=1 log f j (θj ; xij )
2. Estimer θ0,C par

IPM
θC = arg max `(θC , b
b θ1 , . . . , b
θd ; X).
θC ∈ΘC

IPM
θC , b
Les fonctions d’inférence pour les estimateurs sont donc (b θ1 , . . . , b
θd ).

Cette méthode est facile à implémenter puisqu’elle implique typiquement d optimisa-

tions de dimension respective |θj | pour j = 1, . . . , d ainsi qu’une optimisation pour θ0,C .
Les estimés sont prouvables asymptotiquement normaux. Ils sont en revanche sujets à
propagation de l’erreur en cas de mis-spécification des distributions marginales.

Inférence semiparamétrique

Une autre approche, la maximisation de la pseudo-vraisemblance, est suggérée par

Oakes (1994), Genest & Rivest (1995) et adaptée pour les données censurées par Shih &
Louis (1995) consiste à estimer les marginales nonparamétriquement.

Proposition 8.19 (Estimateur du maximum de pseudo-vraisemblance)

À partir d’un échantillon X,

1. Calculer les pseudos-observations basées sur les rangs U

e 1, . . . , U
e n où

n e Rij
U
e ij = Fn,j ( Xij ) =
n+1 n+1

et Rij dénote le rang de Xij parmi X1j , . . . , Xnj .

107
 2. Estimer b
θ0,C par

n n
EMPV
θC
b = arg max ∑ `C (θC ; U e id ) = arg max ∑ log c(U
e i1 , . . . , U e i ; θC ).
θC ∈ΘC i =1 θC ∈ΘC i =1

L’estimateur du maximum de pseudo-vraisemblance n’est pas asymptotiquement effi-

cace, mais est plus robuste et adéquat pour n grand. Cette méthode est également facile
à implémenter, les estimés sont asymptotiquement normaux et crucialement invariants
à des transformations monotones croissantes des données puisque basés sur des rangs.
Il ne faut cependant pas que les marges dépendent de covariables.

La dernière méthode semi-paramétrique est basée sur le méthode des moments. Dans
le cas où θ0,C est unidimensionnel, il est possible d’obtenir pour plusieurs familles de
copules des formules explicites pour le τ de Kendall ou le ρ de Spearman. On peut
alors procéder à l’inversion de $ ou de τ : si la copule est paramétrisée par θ0,C , cela
revient à estimer θbn := τ −1 (τn ). τn est une U-statistique et de ce fait il est possible de
prouver la normalité asymptotique de cette méthode basée sur la méthode δ.

Validation de modèles

Une fois l’estimation complétée, une question cruciale demeure la validation du mo-
dèle. À ce sujet, quelques tests d’adéquations, c’est-à-dire le test de l’hypothèse H0 :
C ∈ Cθ , sont disponibles : voir Genest & al (2009) et la revue de littérature de l’ar-
ticle, de même que Kojadinovic & Yan (2011). Ces tests sont basés sur une technique
d’auto-amorçage (bootstrap) paramétrique ou multiplicatif, et parmi les hypothèses
nécessaires à leur application figurent notamment l’hypothèse d’unicité des rangs (ab-
sence de doublons).
Remarque
L’adéquation ne doit pas être confondue avec la sélection de modèles. Le praticien doit
utiliser un modèle dont les propriétés empiriques ressemblent à celle du modèle choisi.
Peu d’articles se penchent sur la sélection, mais Huard, Évin & Favre (2006) proposent
une alternative bayésienne.

108
 Index
ACP, 36 corrélation, 5
composante principale, 36 Cramér-Wold, 10
test d’isotropie partielle, 39
algorithme espérance-maximisation, 63 dépendance de queue, 99
analyse canonique, 49 dendrogramme, 57
coefficient de corrélation canonique, diagramme
50 P–P, 17
variables canoniques, 50 Q–Q, 17
analyse des correspondances, 53 discrimination
analyse discriminante, voir discrimination fonction de Fisher, 77
approximation linéaire, 76
de Bartlett, 30 logistique, 81
de Box, 29 méthode de Fisher, 78
arbre qualité d’une, 80
classification, 83 dissimilarité, 57
decision, 83
EM
élagage, 83
algorithme, 66
impureté de noeud, 84
données complètes, 65
bornes de Fréchet–Hoeffding, 95 données manquantes, 66
information observée, 70
classification, 73 mélange fini, 69
CEMC, 73 variables latentes, 65
lois multinormales, 76 énergie
cluster, voir partitionnement de données test de normalité, 21
copule, 89
antimonotone, 89 fonction caractéristique, 9
archimédienne, 101 fonction discriminante, 79
Clayton, 104 formule de Lance et Williams, 59
comonotone, 89
Glivenko–Cantelli, 88
densité, 106
elliptique, 99 χ2
empirique, 106 test d’adéquation de, 18
gaussienne, 100 Kolmogorov-Smirnov
Gumbel–Hougaard, 104 test d’adéquation de, 19
indépendance, 89
inférence, 105 loi
Marshall–Olkin, 89 Hotelling, 15
méta, 94 multinormale, 8
de valeurs extrêmes, 104 normale, 8

109
 Wilks, 29 silhouette, 61
Wishart, 13 similarité, 57
Sklar, théorème de, 90
Mahalanobis statistique de Box, 35
distance de, 20, 76
test de normalité, 20 τ de Kendall, 98
MANOVA, 30 théorème de Frobenius, 41
critères, 33 transformation des quantiles, 87
Mardia transformation linéaire, 6
test de normalité, 20 transformée intégrale de probabilité, 87
Marshall-Olkin, 104
matrice variance, 5
de poids, 46 conjointe, 31
de chargement, 45, 46 inter-groupe, 31
de scores, 45, 46 intra-groupe, 31
matrice de variance, 5
Wishart
méthodes hiérarchiques, 58
distribution, voir loi Wishart
barycentre, 58
propriétés, 13
distance complète, 58
distance moyenne, 58
distance unique, 58
médiane, 58
Ward, 59
moyenne, 4
multinomiale, 83

NIPALS, algorithme, 46

partitionnement de données, 57
pseudo-observations, 86

régression
lasso, 44
PLS, 45
ridge, 43
SCAD, 44
sur composantes principales, 45
régression logistique, 82
ρ de Pearson, 96
$ de Spearman, 97

Shapiro–Wilks
test d’adéquation de, 19

110