Calculez une analyse en composantes principales

Un relevé de températures moyennes donne, pour chaque ville, une valeur par saison ou par mois : douze colonnes, et autant de graphiques possibles. Comment voir d’un coup d’œil quelles villes se ressemblent et ce qui les distingue ? L’analyse en composantes principales (ACP, en anglais principal component analysis ou PCA) répond à cette question : elle remplace de nombreuses variables corrélées par quelques variables nouvelles, non corrélées, qui gardent l’essentiel des différences entre individus.

Ce chapitre vous montre comment l’ACP se calcule, sur des exemples assez petits pour être suivis à la main. Le chapitre suivant apprendra à lire le résultat. Notre fil rouge : six villes françaises décrites par leur température moyenne à chaque saison.

L’idée : résumer sans perdre les différences

Commençons par un exemple minuscule. Quatre étudiants ont eu les notes suivantes :

	Français	Philosophie	Maths	Physique
Étudiant 1	16	14	8	10
Étudiant 2	17	15	7	9
Étudiant 3	6	8	18	16
Étudiant 4	7	9	17	15

On aimerait résumer ces quatre notes par deux nombres. Premier essai : la moyenne de français et de maths, et la moyenne de philosophie et de physique. Calculez-les : tous les étudiants obtiennent 12 et 12. Ce résumé est inutile, il efface toutes les différences.

Deuxième essai : la moyenne des matières littéraires (français, philosophie) et celle des matières scientifiques (maths, physique). On obtient 15 et 9 pour l’étudiant 1, 16 et 8 pour le 2, 7 et 17 pour le 3, 8 et 16 pour le 4. Cette fois, deux profils sautent aux yeux : deux « littéraires » et deux « scientifiques ».

Deux résumés en deux nombres des mêmes notes : le premier écrase tout le monde au même point, le second fait apparaître les deux profils. L'ACP cherche le second type de résumé.

Le bon résumé est donc celui qui conserve le mieux la variabilité entre les individus. L’ACP trouve automatiquement ces bonnes combinaisons, avec deux objectifs :

1. transformer un grand nombre de variables corrélées en un petit nombre de variables non corrélées, les composantes principales ;
2. choisir celles qui montrent le mieux la variabilité du nuage de points.

Pensez à une photographie d’un objet en trois dimensions. Selon l’angle, on reconnaît l’objet ou pas : une tasse vue de dessus n’est qu’un cercle, vue de profil on voit l’anse. L’ACP cherche l’angle de vue qui déforme le moins le nuage de points.

Deux nuages et une préparation

Les données de l’ACP forment un tableau de $I$ individus (les lignes) décrits par $K$ variables quantitatives (les colonnes). Ce tableau se lit de deux façons :

• le nuage des individus : chaque individu est un point de $\mathbb{R}^K$ ; on cherche quels individus se ressemblent ;
• le nuage des variables : chaque variable est un point de $\mathbb{R}^I$ ; on cherche quelles variables sont liées.

Avant tout calcul, on centre toujours les données : on retire à chaque variable sa moyenne, $x_{ik} \leftarrow x_{ik} - \bar x_k$, où $x_{ik}$ est la valeur de la variable $k$ pour l’individu $i$ et $\bar x_k$ la moyenne de la variable $k$. L’origine se place alors au centre de gravité du nuage, sans en changer la forme.

On réduit aussi, le plus souvent : $x_{ik} \leftarrow (x_{ik} - \bar x_k)/s_k$, où $s_k$ est l’écart-type de la variable $k$. C’est indispensable dès que les variables n’ont pas la même unité ou le même ordre de grandeur : sans réduction, une variable exprimée en euros (variance énorme) écraserait une variable exprimée en années. On parle alors d’ACP normée. Dans la suite, on note $Z$ le tableau centré-réduit.

Attention

L’ACP ne s’applique qu’à des variables quantitatives. Pour des variables qualitatives, il existe des méthodes cousines, l’AFC et l’ACM, présentées au chapitre suivant.

Le critère : l’inertie projetée

Pour mesurer la dispersion d’un nuage à plusieurs dimensions, on généralise la variance : c’est l’inertie, la moyenne des carrés des distances des points au centre de gravité $G$.

$$I_{\text{tot}} = \frac{1}{I}\sum_{i=1}^{I} d^2(M_i, G) = \sum_{k=1}^{K} \operatorname{Var}(x_k)$$

où $M_i$ est le point représentant l’individu $i$. L’inertie totale est donc la somme des variances des variables. En ACP normée, chaque variable a une variance de 1, et l’inertie totale vaut simplement $K$.

Projetons maintenant chaque point $M_i$ sur un axe passant par l’origine $O$ (le centre de gravité, puisque les données sont centrées). Appelons $H_i$ son projeté.

Le triangle $O H_i M_i$ est rectangle en $H_i$, donc par Pythagore : $OM_i^2 = OH_i^2 + M_iH_i^2$. La longueur $OM_i$ ne dépend pas de l’axe choisi. Par conséquent, minimiser la déformation $\sum_i M_iH_i^2$ revient exactement à maximiser l’inertie projetée $\sum_i OH_i^2$ : l’axe qui déforme le moins le nuage est celui sur lequel les projections sont les plus étalées.

L’ACP procède alors pas à pas :

1. trouver l’axe $u_1$ d’inertie projetée maximale ;
2. trouver l’axe $u_2$, orthogonal à $u_1$, d’inertie projetée maximale ; le plan $(u_1, u_2)$ est le meilleur plan ;
3. continuer avec $u_3$ orthogonal aux deux premiers, et ainsi de suite jusqu’à $K$ axes.

Le calcul : valeurs et vecteurs propres

Rappelons qu’un vecteur non nul $v$ est un vecteur propre d’une matrice carrée $A$ si $Av = \lambda v$ : la matrice ne fait que l’allonger ou le raccourcir, d’un facteur $\lambda$ appelé valeur propre. On montre que les axes cherchés sont les vecteurs propres de la matrice des corrélations. Le calcul tient en quatre étapes.

1. Calculer la matrice des corrélations $R = \frac{1}{I} Z^\top Z$, où $Z$ est le tableau centré-réduit. C’est une matrice $K \times K$ symétrique, dont le terme $(k, l)$ est la corrélation entre les variables $k$ et $l$.
2. Calculer ses valeurs propres $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_K \geq 0$ et ses vecteurs propres unitaires $u_1, \dots, u_K$. Comme $R$ est symétrique, ces vecteurs sont orthogonaux deux à deux.
3. L’axe $s$ est dirigé par $u_s$. La composante principale $F_s = Z u_s$ donne la coordonnée de chaque individu sur cet axe. Elle est de moyenne nulle et de variance $\lambda_s$.
4. L’inertie expliquée par l’axe $s$ vaut $\lambda_s$ ; celle du plan $(1, 2)$ vaut $\lambda_1 + \lambda_2$. La somme de toutes les valeurs propres est l’inertie totale, $K$ en ACP normée.

Exemple

Prenons deux variables standardisées de corrélation $r = 0{,}8$. La matrice des corrélations est

$$R = \begin{pmatrix} 1 & 0{,}8 \\ 0{,}8 & 1 \end{pmatrix}$$

Essayons $u = \tfrac{1}{\sqrt 2}(1, 1)$ : $R u = \tfrac{1}{\sqrt 2}(1{,}8 ;\, 1{,}8) = 1{,}8\, u$. C’est un vecteur propre, de valeur propre $\lambda_1 = 1{,}8$. De même, $\tfrac{1}{\sqrt 2}(1, -1)$ est vecteur propre de valeur propre $\lambda_2 = 0{,}2$. On vérifie que $\lambda_1 + \lambda_2 = 2 = K$.

Le premier axe explique $1{,}8 / 2 = 90\,\%$ de l’inertie. La première composante, $F_1 = \tfrac{1}{\sqrt 2}(z_1 + z_2)$, est une moyenne (à un facteur près) des deux variables : quand deux variables sont très corrélées, l’ACP les fusionne en une seule.

Plus généralement, pour deux variables de corrélation $r$ positive, les valeurs propres valent $1 + r$ et $1 - r$ : plus la corrélation est forte, plus le premier axe résume à lui seul.

Nuage simulé avec r = 0,8. Le premier axe suit l'allongement du nuage et porte 90 % de l'inertie ; la projection Hᵢ de chaque point sur cet axe est sa première composante principale, au facteur près.

Et nos quatre étudiants ? Leurs notes sont presque parfaitement corrélées (plus de 0,99 en valeur absolue). L’ACP normée donne $\lambda_1 \approx 3{,}996$, soit 99,9 % de l’inertie sur le premier axe, dirigé par $u_1 = (0{,}5 ;\, 0{,}5 ;\, -0{,}5 ;\, -0{,}5)$. La première composante vaut donc la moitié de (français + philosophie − maths − physique), en notes standardisées : c’est exactement l’opposition « littéraire contre scientifique » que nous avions devinée.

Le fil rouge : six villes, quatre saisons

Voici maintenant un exemple plus réaliste. Les températures moyennes (en °C, valeurs arrondies construites pour l’exemple) de six villes sont :

Ville	Hiver	Printemps	Été	Automne
Brest	7,0	10,5	16,5	12,0
Rennes	5,5	11,0	18,0	12,0
Strasbourg	1,5	10,5	19,5	10,0
Lyon	3,5	12,0	21,5	12,5
Nice	9,5	13,5	23,0	17,0
Marseille	8,0	14,0	24,0	16,0

Les lignes montent toutes vers l'été : c'est la chaleur générale. Mais elles se croisent : Strasbourg est la plus froide en hiver et plus chaude que Brest en été. Ce croisement est la seconde source de différences.

Toutes les variables sont en degrés, mais leurs variances diffèrent ; on fait une ACP normée. La matrice des corrélations vaut (arrondie au millième) :

	Hiver	Printemps	Été	Automne
Hiver	1	0,658	0,345	0,870
Printemps	0,658	1	0,926	0,933
Été	0,345	0,926	1	0,757
Automne	0,870	0,933	0,757	1

Toutes les corrélations sont positives : une ville chaude à une saison a tendance à l’être aux autres. Mais l’hiver et l’été ne sont liés qu’à 0,345, signe qu’il existe une seconde source de différences. Les valeurs propres de cette matrice sont :

Axe	Valeur propre	% d’inertie	% cumulé
1	3,273	81,8 %	81,8 %
2	0,714	17,9 %	99,7 %
3	0,011	0,3 %	99,9 %
4	0,002	0,1 %	100 %

Les vecteurs propres des deux premiers axes sont $u_1 \approx (0{,}437 ;\, 0{,}541 ;\, 0{,}468 ;\, 0{,}545)$ et $u_2 \approx (0{,}725 ;\, -0{,}222 ;\, -0{,}627 ;\, 0{,}178)$, dans l’ordre hiver, printemps, été, automne. Le premier axe donne un poids positif et comparable aux quatre saisons : il mesure si une ville est chaude toute l’année. Le second oppose l’hiver (poids positif) à l’été (poids négatif) : il mesure la douceur de l’hiver par rapport à l’été, autrement dit l’amplitude thermique. Le chapitre suivant confirmera cette lecture.

u₁ pèse à peu près autant sur les quatre saisons (chaleur générale) ; u₂ oppose l'hiver (poids positif) à l'été (poids négatif).

Combien d’axes garder ?

Trois critères coexistent :

• le pourcentage d’inertie cumulé : $\tau_s = \lambda_s / \sum_k \lambda_k$ ; les axes étant orthogonaux, les pourcentages s’additionnent. On garde assez d’axes pour atteindre un seuil fixé selon l’usage (souvent 80 à 90 %) ;
• le coude dans le diagramme des valeurs propres (l’éboulis) : on s’arrête là où la décroissance devient brutalement plus lente ;
• le critère de Kaiser, en ACP normée : on garde les axes de valeur propre supérieure à 1, c’est-à-dire ceux qui résument plus d’information qu’une variable d’origine à elle seule.

Pour nos villes, le plan $(1, 2)$ explique 99,7 % de l’inertie : il est inutile d’aller plus loin. L’exemple classique de 15 grandes villes françaises décrites par leurs 12 températures mensuelles donne un résultat très proche : près de 80 % de l’inertie sur le premier axe, près de 19 % sur le second, soit presque 99 % pour le plan.

Attention

Les critères ne sont pas toujours d’accord. Ici, Kaiser ne garderait que le premier axe ($\lambda_2 = 0{,}714$ est inférieur à 1), alors que le second porte 17,9 % de l’inertie et une information très lisible, l’amplitude thermique. Avec peu de variables, Kaiser est trop sévère. Aucun critère ne remplace l’interprétation.

La courbe s'aplatit après l'axe 2 : le coude désigne deux axes (99,7 % cumulés), alors que Kaiser n'en garderait qu'un : la seconde valeur propre, 0,714, passe juste sous la ligne λ = 1.

À retenir

L’ACP normée diagonalise la matrice des corrélations. Les vecteurs propres donnent les axes, les valeurs propres l’inertie expliquée par chacun, et les composantes principales $F_s = Z u_s$ les coordonnées des individus. On garde les premiers axes selon le pourcentage d’inertie cumulé, le coude et, avec prudence, le critère de Kaiser.

En pratique : la SVD et scikit-learn

Les logiciels ne diagonalisent pas $R$ directement : ils passent par la décomposition en valeurs singulières (SVD) du tableau centré-réduit, $Z = U \Sigma V^\top$, où $U$ et $V$ ont des colonnes orthonormées et $\Sigma$ est diagonale, de termes $\sigma_1 \geq \sigma_2 \geq \dots$ (les valeurs singulières). On en déduit $Z^\top Z = V \Sigma^2 V^\top$ : les axes principaux sont les colonnes de $V$, et les valeurs propres valent $\lambda_s = \sigma_s^2 / I$. Ce calcul est plus stable numériquement, et il évite de former la matrice $Z^\top Z$.

import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

villes = pd.DataFrame(
    {"hiver": [7.0, 5.5, 1.5, 3.5, 9.5, 8.0], "printemps": [10.5, 11.0, 10.5, 12.0, 13.5, 14.0],
     "ete": [16.5, 18.0, 19.5, 21.5, 23.0, 24.0], "automne": [12.0, 12.0, 10.0, 12.5, 17.0, 16.0]},
    index=["Brest", "Rennes", "Strasbourg", "Lyon", "Nice", "Marseille"])

Z = StandardScaler().fit_transform(villes)     # centrer-réduire (écart-type divisé par n)
acp = PCA().fit(Z)
print(acp.explained_variance_ratio_.round(3))  # [0.818 0.179 0.003 0.001]
print(acp.components_[:2].round(3))            # u1 et u2 (au signe près)
F = acp.transform(Z)                           # coordonnées des villes sur les axes

Piège

Deux détails de scikit-learn surprennent souvent. D’abord, explained_variance_ divise par $I - 1$ et non par $I$ : sur 6 villes, il affiche $3{,}273 \times 6/5 \approx 3{,}93$ au lieu de 3,273. Les pourcentages (explained_variance_ratio_) ne sont pas affectés. Ensuite, le signe d’un vecteur propre est arbitraire : $u_1$ et $-u_1$ définissent le même axe. Un logiciel peut donc afficher un axe « à l’envers » ; l’interprétation reste la même, il suffit de retourner la lecture.

En résumé

• L’ACP remplace $K$ variables quantitatives corrélées par des composantes principales non corrélées qui conservent au mieux la variabilité du nuage.
• On centre toujours les données, et on les réduit dès que les unités ou les ordres de grandeur diffèrent (ACP normée).
• Le critère est l’inertie projetée : par Pythagore, la maximiser revient à déformer le nuage le moins possible.
• Les axes sont les vecteurs propres de la matrice des corrélations ; chaque valeur propre $\lambda_s$ est l’inertie (la variance) expliquée par l’axe $s$ ; leur somme vaut $K$ en ACP normée.
• Sur les six villes, le premier axe (81,8 %) mesure la chaleur générale, le second (17,9 %) l’amplitude thermique.
• On choisit le nombre d’axes par le pourcentage cumulé, le coude ou le critère de Kaiser ($\lambda > 1$), qui peut être trop sévère.

Exercices

Exercice · Reconnaître un vecteur propre

facile

Soit $A = \begin{pmatrix} 2 & 3 \\ 2 & 1 \end{pmatrix}$. Les vecteurs $v = (3, 2)$ et $w = (2, 3)$ sont-ils des vecteurs propres de $A$ ? Si oui, avec quelle valeur propre ? Le vecteur $(1, -1)$ est aussi un vecteur propre de $A$ : est-il orthogonal à $v$ ? Qu’en concluez-vous ?

Voir la solution

• $A v = (2 \times 3 + 3 \times 2 ;\, 2 \times 3 + 1 \times 2) = (12, 8) = 4\,(3, 2)$ : $v$ est vecteur propre, de valeur propre 4.
• $A w = (2 \times 2 + 3 \times 3 ;\, 2 \times 2 + 1 \times 3) = (13, 7)$, qui n’est pas proportionnel à $(2, 3)$ : $w$ n’est pas vecteur propre.
• $A(1, -1) = (2 - 3 ;\, 2 - 1) = (-1, 1) = -1 \times (1, -1)$ : vecteur propre de valeur propre $-1$. Son produit scalaire avec $v$ vaut $3 \times 1 + 2 \times (-1) = 1 \neq 0$ : ils ne sont pas orthogonaux.

Conclusion : les vecteurs propres d’une matrice quelconque ne sont pas forcément orthogonaux. Ils le sont pour une matrice symétrique, comme la matrice des corrélations de l’ACP : c’est ce qui garantit des axes perpendiculaires.

En code

Python

import numpy as np

A = np.array([[2, 3], [2, 1]])
v, w, t = np.array([3, 2]), np.array([2, 3]), np.array([1, -1])

print(A @ v, A @ w, A @ t)     # [12 8] = 4 v ; [13 7] non proportionnel ; [-1 1] = -1 t
print(v @ t)                   # 1 : v et t ne sont pas orthogonaux
print(np.linalg.eigvals(A).real)   # [ 4. -1.]

# Une matrice symétrique, elle, a des vecteurs propres orthogonaux
S = np.array([[1, 0.8], [0.8, 1]])
_, vecteurs = np.linalg.eigh(S)
print(round(vecteurs[:, 0] @ vecteurs[:, 1], 10))   # 0.0

R

A <- matrix(c(2, 3,
              2, 1), nrow = 2, byrow = TRUE)
v <- c(3, 2); w <- c(2, 3); t1 <- c(1, -1)

A %*% v; A %*% w; A %*% t1   # (12, 8) = 4 v ; (13, 7) non proportionnel ; (-1, 1) = -1 t1
sum(v * t1)                  # 1 : v et t1 ne sont pas orthogonaux
eigen(A)$values              # 4 -1

# Une matrice symétrique, elle, a des vecteurs propres orthogonaux
S <- matrix(c(1, 0.8, 0.8, 1), 2)
vecteurs <- eigen(S)$vectors
round(sum(vecteurs[, 1] * vecteurs[, 2]), 10)   # 0

Exercice · Une ACP à deux variables

moyen

Deux variables standardisées ont une corrélation $r = -0{,}6$.

1. Écrivez la matrice des corrélations et vérifiez que $\tfrac{1}{\sqrt 2}(1, -1)$ et $\tfrac{1}{\sqrt 2}(1, 1)$ sont des vecteurs propres. Quelles sont les valeurs propres ?
2. Quel pourcentage d’inertie le premier axe explique-t-il ?
3. Un individu a pour valeurs standardisées $z = (1, -1)$. Calculez ses coordonnées $F_1$ et $F_2$.

Voir la solution

1. $R = \begin{pmatrix} 1 & -0{,}6 \\ -0{,}6 & 1 \end{pmatrix}$. Pour $a = \tfrac{1}{\sqrt 2}(1, -1)$ : $R a = \tfrac{1}{\sqrt 2}(1 + 0{,}6 ;\, -0{,}6 - 1) = 1{,}6\, a$. Pour $b = \tfrac{1}{\sqrt 2}(1, 1)$ : $R b = \tfrac{1}{\sqrt 2}(0{,}4 ;\, 0{,}4) = 0{,}4\, b$. Les valeurs propres sont $\lambda_1 = 1{,}6$ (axe $u_1 = a$) et $\lambda_2 = 0{,}4$ (axe $u_2 = b$). Leur somme vaut bien 2.
2. $1{,}6 / 2 = 80\,\%$. Avec une corrélation négative, le premier axe oppose les deux variables au lieu de les additionner.
3. $F_1 = z \cdot u_1 = \tfrac{1}{\sqrt 2}(1 \times 1 + (-1) \times (-1)) = \tfrac{2}{\sqrt 2} = \sqrt 2 \approx 1{,}414$ et $F_2 = z \cdot u_2 = \tfrac{1}{\sqrt 2}(1 - 1) = 0$. L’individu est entièrement porté par le premier axe : il est fort sur la première variable et faible sur la seconde, exactement l’opposition que résume $u_1$.

En code

Python

import numpy as np

R = np.array([[1, -0.6], [-0.6, 1]])
a = np.array([1, -1]) / np.sqrt(2)
b = np.array([1, 1]) / np.sqrt(2)
print(R @ a / a, R @ b / b)          # [1.6 1.6]  [0.4 0.4] : valeurs propres 1.6 et 0.4

lambdas = np.array([1.6, 0.4])
print(lambdas[0] / lambdas.sum())    # 0.8 : 80 % de l'inertie sur l'axe 1

z = np.array([1, -1])
print(round(z @ a, 3), z @ b)        # F1 = 1.414, F2 = 0.0

R

R <- matrix(c(1, -0.6, -0.6, 1), 2)
a <- c(1, -1) / sqrt(2)
b <- c(1, 1) / sqrt(2)
drop(R %*% a) / a; drop(R %*% b) / b   # 1.6 1.6 ; 0.4 0.4 : valeurs propres

lambdas <- eigen(R)$values             # 1.6 0.4
lambdas[1] / sum(lambdas)              # 0.8 : 80 % de l'inertie sur l'axe 1

z <- c(1, -1)
round(c(F1 = sum(z * a), F2 = sum(z * b)), 3)   # F1 = 1.414, F2 = 0

Exercice · Pourquoi la variance d'une composante vaut λ

difficile

Soit $Z$ un tableau centré-réduit ($I$ lignes, $K$ colonnes), $R = \frac{1}{I} Z^\top Z$, et $u_s$ un vecteur propre unitaire de $R$ de valeur propre $\lambda_s$. On pose $F_s = Z u_s$.

1. Montrez que $F_s$ est de moyenne nulle.
2. Montrez que sa variance vaut $\lambda_s$.
3. Déduisez-en que la somme des valeurs propres vaut $K$.

Voir la solution

1. La moyenne de $F_s$ est $\frac{1}{I}\sum_i F_s(i) = \frac{1}{I}\sum_i \sum_k z_{ik} u_{ks} = \sum_k u_{ks} \left(\frac{1}{I}\sum_i z_{ik}\right)$. Chaque colonne de $Z$ est centrée, donc chaque parenthèse est nulle : la moyenne vaut 0.
2. Comme la moyenne est nulle, la variance est la moyenne des carrés : $\operatorname{Var}(F_s) = \frac{1}{I} F_s^\top F_s = \frac{1}{I} u_s^\top Z^\top Z u_s = u_s^\top R u_s$. Or $R u_s = \lambda_s u_s$, donc $u_s^\top R u_s = \lambda_s u_s^\top u_s = \lambda_s$, puisque $u_s$ est unitaire.
3. La somme des valeurs propres d’une matrice est égale à sa trace (la somme de ses termes diagonaux). La diagonale de $R$ contient les corrélations de chaque variable avec elle-même, toutes égales à 1 : la trace vaut $K$. C’est bien l’inertie totale d’une ACP normée, répartie entre les axes.

En code

Python

import numpy as np

# Vérification sur les six villes du chapitre
X = np.array([[7.0, 10.5, 16.5, 12.0], [5.5, 11.0, 18.0, 12.0], [1.5, 10.5, 19.5, 10.0],
              [3.5, 12.0, 21.5, 12.5], [9.5, 13.5, 23.0, 17.0], [8.0, 14.0, 24.0, 16.0]])
I, K = X.shape
Z = (X - X.mean(axis=0)) / X.std(axis=0)       # centré-réduit (écart-type divisé par I)
R = Z.T @ Z / I
lambdas, U = np.linalg.eigh(R)                 # vecteurs propres unitaires en colonnes
F = Z @ U                                      # toutes les composantes d'un coup

print(F.mean(axis=0).round(10))                # [0. -0. -0. -0.] : moyennes nulles
print(F.var(axis=0).round(3))                  # [0.002 0.011 0.714 3.273]
print(lambdas.round(3))                        # les mêmes : Var(F_s) = lambda_s
print(round(lambdas.sum(), 6), round(np.trace(R), 6))   # 4.0  4.0 : somme = trace = K

R

# Vérification sur les six villes du chapitre
X <- matrix(c(7.0, 10.5, 16.5, 12.0,  5.5, 11.0, 18.0, 12.0,  1.5, 10.5, 19.5, 10.0,
              3.5, 12.0, 21.5, 12.5,  9.5, 13.5, 23.0, 17.0,  8.0, 14.0, 24.0, 16.0),
            ncol = 4, byrow = TRUE)
I <- nrow(X); K <- ncol(X)
ecart_n <- function(x) sqrt(mean((x - mean(x))^2))            # écart-type divisé par I
Z <- scale(X, center = TRUE, scale = apply(X, 2, ecart_n))
R <- t(Z) %*% Z / I
e <- eigen(R)                                  # vecteurs propres unitaires en colonnes
F <- Z %*% e$vectors                           # toutes les composantes d'un coup

round(colMeans(F), 10)                         # 0 0 0 0 : moyennes nulles
round(colMeans(F^2), 3)                        # variances : 3.273 0.714 0.011 0.002
round(e$values, 3)                             # les mêmes : Var(F_s) = lambda_s
c(sum(e$values), sum(diag(R)))                 # 4 4 : somme = trace = K