Analyse en composantes principales¶
Voyons un exemple avec des données simulées.
%matplotlib nbagg
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch
On simule deux jeux de données de même taille. Le premier tableau de données est composé de 20 réalisations i.i.d. de la loi normale $\mathcal N(0,0.3)$, le deuxième contient 20 réalisations i.i.d. de la loi normale $\mathcal N(0.5,.2)$.
np.random.seed(234)
jeu1 = pd.DataFrame(np.sqrt(.3)*np.random.randn(20,3),columns=list('ABC'))
print 'jeu1 :', jeu1.shape, '\n', jeu1.head(), '\n'
jeu2 = pd.DataFrame(np.sqrt(.2)*np.random.randn(20,3)+.5,columns=list('ABC'))
print 'jeu2 :', jeu2.shape, '\n', jeu2.head()
print jeu1.describe()
print jeu2.describe()
Voici quelques représentations graphiques des données. D'abord variable par variable :
plt.figure(figsize=(8,8))
plt.subplot(311)
plt.boxplot([jeu1['A'],jeu2['A']])
plt.ylabel('A')
plt.subplot(312)
plt.boxplot([jeu1['B'],jeu2['B']])
plt.ylabel('B')
plt.subplot(313)
plt.boxplot([jeu1['C'],jeu2['C']])
plt.ylabel('C')
plt.figure(figsize=(8,8))
plt.subplot(311)
plt.plot(jeu1['A'], np.zeros(20), 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['A'], np.zeros(20), '^', markersize=8, color='red', label='jeu2')
plt.legend(loc=0)
plt.subplot(312)
plt.plot(jeu1['B'], np.zeros(20), 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['B'], np.zeros(20), '^', markersize=8, color='red', label='jeu2')
plt.subplot(313)
plt.plot(jeu1['C'], np.zeros(20), 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['C'], np.zeros(20), '^', markersize=8, color='red', label='jeu2')
Pour les trois variables on observe beaucoup de chevauchements des deux groupes. Autrement dit, on ne peut pas utiliser une des variables pour séparer les deux jeux de données.
Voyons avec les nuages des points en dimension 2 :
plt.figure(figsize=(8,8))
plt.subplot(311)
plt.plot(jeu1['A'], jeu1['B'], 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['A'], jeu2['B'], '^', markersize=8, color='red', label='jeu2')
plt.title('Nuage des points de jeu1 et jeu2')
plt.legend(loc=0)
plt.subplot(312)
plt.plot(jeu1['A'], jeu1['C'], 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['A'], jeu2['C'], '^', markersize=8, color='red', label='jeu2')
plt.subplot(313)
plt.plot(jeu1['B'], jeu1['C'], 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['B'], jeu2['C'], '^', markersize=8, color='red', label='jeu2')
Les deux jeux de données sont légèrement mieux séparés qu'en dimension 1, mais la séparation est loin d'être optimale.
Enfin, traçons le nuage des points en 3D :
fig = plt.figure()
fig.add_subplot(111, projection='3d')
plt.plot(jeu1['A'], jeu1['B'], jeu1['C'], 'o', markersize=8, color='blue', label='jeu1')
plt.plot(jeu2['A'], jeu2['B'], jeu2['C'], '^', markersize=8, color='red', label='jeu2')
plt.title('Nuage des points de jeu1 et jeu2')
plt.legend(loc=0)
plt.xlabel('A')
plt.ylabel('B')
Noter que la figure ci-dessus est interactive : vous pouvez la faire tourner en cliquant et bougeant la souris. Essayez-le !
On dirait que les points rouges se trouvent plutôt en haut à droite, alors que les points bleus se trouvent en bas à gauche. Il semble qu'on pourrait introduire un plan qui sépare assez bien les deux jeux de données. Mais comment trouver ce plan ? C'est par l'analyse en composantes principales (ACP) !
En fait, l'ACP est une technique de réduction de dimension des données. Plus précisément, on cherche un sous-espace vectoriel de petite dimension qui explique au mieux les données, c'est-à-dire qui préserve un maximum de la variance des données quand on les projette dans cet sous-espace vectoriel.
L'ACP repose sur une SVD des données complètes. On crée alors un seul tableau de données :
df = pd.concat([jeu1,jeu2],ignore_index=True)
print 'df :', df.shape, '\n', df.head()
df.describe()
Rappelons que la SVD d'une matrice $X$ de taille $n\times p$ consiste à calculer les valeurs singulières $\lambda=(\lambda_1,\dots,\lambda_p)$ (triées par ordre décroissant) et les matrices unitaires $U=[u_1,\dots,u_n]$ et $V=[v_1,\dots,v_p]$, tels que $$X=\sum_{k=1}^p \lambda_ku_kv_k^T = U\text{diag}(\lambda)V^T.$$ Les colonnes $v_k$ de $V$ sont les axes principaux et les vecteurs $\lambda_k u_k$ sont les composantes principales.
En Python, c'est la fonction np.linalg.svd qui effectue la SVD d'une matrice. Elle renvoie les valeurs singulières $\lambda$ et les matrices unitaires $U$ et $V$.
Pour l'ACP, nous effectuons la SVD du tableau des données centré.
moy = df.mean()
df_centre = df
df_centre['A'] =df_centre['A'] - moy[0]
df_centre['B'] =df_centre['B'] - moy[1]
df_centre['C'] =df_centre['C'] - moy[2]
U, lamb ,V = np.linalg.svd(df_centre)
print 'U : dimensions :', U.shape, '\n' , U, '\n lambda :\n', lamb, '\n V :dimensions :', V.shape, '\n', V, '\n'
Afin d'illustrer les axes principaux, on les ajoute au graphique du nuage des points. Por cela, nous définissons d'abord une classe pour tracer des flêches :
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
self._verts3d = xs, ys, zs
def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
FancyArrowPatch.draw(self, renderer)
Voici le nuage des points avec les directions des axes principaux :
fig2 = plt.figure(figsize=(7,7))
ax = fig2.add_subplot(111, projection='3d')
plt.plot(df['A'],df['B'],df['C'], 'o', markersize=8, color='green',alpha=0.7)
plt.plot([moy[0]], [moy[1]], [moy[2]], 'o', markersize=10, color='red')
for i in range(V.shape[0]):
print V[i,0], V[i,1], V[i,2]
a = Arrow3D([moy[0], V[i,0]], [moy[1], V[i,1]], [moy[2], V[i,2]], mutation_scale=20, lw=3, arrowstyle="-|>", color="r")
ax.add_artist(a)
plt.title('Eigenvectors')
plt.xlabel('A')
plt.ylabel('B')
Calculons les composantes principales :
df_proj = pd.DataFrame(U[:,[0,1,2]], columns=['CP1','CP2','CP3'] )
df_proj['CP1'] = lamb[0]*df_proj['CP1']
df_proj['CP2'] = lamb[1]*df_proj['CP2']
df_proj['CP3'] = lamb[2]*df_proj['CP3']
df_proj.head()
Représentons la première composante principale et comparez aux graphiques obtenus plus haut :
plt.figure()
plt.subplot(211)
plt.boxplot([df_proj.ix[0:19,'CP1'],df_proj.ix[20:39,'CP1']])
plt.ylabel('1er axe principal')
plt.subplot(212)
plt.plot(df_proj.ix[0:19,'CP1'],np.zeros(20), 'o', markersize=8, color='blue', label='jeu2')
plt.plot(df_proj.ix[20:39,'CP1'],np.zeros(20), '^', markersize=8, color='red', label='jeu2')
plt.xlabel('1er axe principal')
Ensuite, représentons les données en dimension 2, c'est-à-dire par le nuage des points défini par les deux premières composantes prinicpales :
plt.figure()
plt.plot(df_proj.ix[0:19,'CP1'],df_proj.ix[0:19,'CP2'], 'o', markersize=8, color='blue', label='jeu2')
plt.plot(df_proj.ix[20:39,'CP1'],df_proj.ix[20:39,'CP2'], '^', markersize=8, color='red', label='jeu2')
plt.xlabel('1er axe principal')
plt.ylabel('2e axe principal')
plt.legend()
plt.title('ACP avec 2 composantes principales')
Comparé aux nuages de points des variables d'origines, les deux jeu de données sont beaucoup mieux séparés. Ceci est dû au fait que l'ACP cherche à préserver un maximum de la variance des données.
Calculons la variance ou l'inertie totale dans les données définie par $$\mathcal I = \frac1n\sum_{i=1}^n\|X_i\|^2$$ où les vecteurs $X_i$ sont supposés centrés.
n = df_centre.shape[0]
inertie_tot = ((df_centre**2).sum()).sum()/n
print 'inertie totale : ', inertie_tot
Dans le cours, on a montré qu'on peut calculer l'inertie également par la formule suivante qui fait intervenir les valeurs singulières : $$\mathcal I= \frac1n\sum_{k=1}^p\lambda^2_k$$
Vérifions que cette formule donne le même résultat :
sum(lamb**2)/n
L'inertie du nuage projeté sur le premier axe principal est donné par $$\mathcal I_1=\frac{\lambda_1^2}{n}$$ et la proportion de l'inertie (ou de la variance) expliquée par la première composante prinicpale est donnée par $$\frac{\lambda_1^2}{\sum_{k=1}^p\lambda^2_k}$$ Ici on trouve pour $I_1$ :
lamb[0]**2/n
et pour la proportion de la variance expliquée on a :
lamb[0]**2/sum(lamb**2)
ça veut dire que pour nos données la première composante principale explique presque la moitié de la variance de nos données d'origine. Pour comparaison, la part de variabilité expliquée par chacune des variables d'origine est entre 0.25 et 0.42 :
df.var()*(n-1)/sum(lamb**2)
Inertie du nuage projeté sur le plan engendré par les deux premiers axes principaux est donnée par $$\mathcal I_2 =\frac{\lambda_1^2+\lambda_2^2}{n}$$ et la proportion de la variance (ou inertie) expliquée par les deux premières composantes prinicpales est donnée par $$\frac{\lambda_1^2+\lambda_2^2}{\sum_{k=1}^p\lambda^2_k}$$
Ici on trouve pour $I_2$ :
sum(lamb[0:2]**2)/n
Et pour la part de la variance expliquée par les deux premières composantes principales :
sum(lamb[0:2]**2)/sum(lamb**2)
Le module sklearn contient une fonction pour faire l'ACP qui est plus commode à utiliser :
from sklearn.decomposition import PCA
sklearn_pca = PCA()
df_transf = sklearn_pca.fit_transform(df) # pas nécessaire de centrer les colonnes ; renvoie les comp. princ.
df_transf
plt.figure()
plt.plot(df_transf[0:19,0],df_transf[0:19,1], 'o', markersize=7, color='blue', label='jeu1')
plt.plot(df_transf[20:39,0], df_transf[20:39,1], '^', markersize=7, color='red', label='jeu2')
plt.xlabel('1er axe principal')
plt.ylabel('2e axe principal')
plt.legend()
plt.title('ACP avec 2 composantes principales')
On accède facilement à la variance (ou inertie) expliquée par chacune des composantes principales, et ces valeurs premettent de retrouver les valeurs singulières. On vérifie qu'on retrouve exactement les mêmes valeurs qu'avec l'autre méthode :
inerties = sklearn_pca.explained_variance_ # inertie expliquée par chaque composante principale
print "Inertie d'une, deux, trois composantes principales :\n", np.cumsum(inerties)
print 'valeurs singulieres :', np.sqrt(inerties*n)
Exercice - Données Iris¶
Un des jeux de données les plus célèbres en statistique est le jeu de données Iris. Ces données quantifient les variations de morphologie des fleurs d'Iris de trois espèces, et elles étaient présentées pour la première fois en 1936 par le statisticien Ronald Fisher.
Un peu de culture générale : Ronald Fisher (1890-1962) est un biologiste et statisticien britannique. Il est considéré par Richard Dawkins comme le plus grand des successeurs de Darwin et par Anders Hald comme l'homme qui a - presque à lui seul - fondé les statistiques modernes. Bradley Efron le considère comme le statisticien le plus important du XXe siècle. Dans le domaine des statistiques, il a introduit de nombreux concepts clés tels que le maximum de vraisemblance et l'information de Fisher.
Le jeu de données Iris comprend 50 échantillons de chacune des trois espèces de Iris (Iris setosa, Iris virginica et Iris versicolor). Quatre caractéristiques ont été mesurées à partir de chaque échantillon : la longueur et la largeur des sépales et les pétales, en centimètres.
Le jeu de données est disponible à l'adresse http://www.proba.jussieu.fr/pageperso/rebafka/irisdata.csv
- Importez les données.
- Familiarisez-vous avec les données. Vérifiez qu'il y a bien le même nombre d'observations pour chaque espèce.
- Tracer des nuages des points (utilisez des couleurs différentes pour les différentes espèces) et des boxplots pour comparer les différentes espèces.
- Effectuez une ACP et représentez graphiquement les résultats.
- Calculez l'inertie totale et l'inertie des données projetées.