# Analyse en composantes principales¶

Voyons un exemple avec des données simulées.

In [1]:
%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

/Users/tabea/anaconda/lib/python2.7/site-packages/IPython/kernel/__init__.py:13: ShimWarning: The IPython.kernel package has been deprecated. You should import from ipykernel or jupyter_client instead.
"You should import from ipykernel or jupyter_client instead.", ShimWarning)


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)$.

In [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()

jeu1 : (20, 3)
A         B         C
0  0.448471 -0.571576  0.192196
1  0.504769 -0.047861 -1.713761
2 -0.531144  0.511938  0.024027
3  0.780623 -0.305116  0.507643
4 -0.703031  0.600445 -1.058459

jeu2 : (20, 3)
A         B         C
0  0.216355  1.095396  0.156774
1  0.903916  1.373958 -0.132441
2  0.329976  0.570340  0.736356
3  0.781577 -0.122020  0.762135
4  1.142091  0.946711  0.794450

In [3]:
print jeu1.describe()
print jeu2.describe()

               A          B          C
count  20.000000  20.000000  20.000000
mean    0.204432   0.063061  -0.081521
std     0.476635   0.523045   0.587254
min    -0.703031  -0.626166  -1.713761
25%    -0.098464  -0.355382  -0.171024
50%     0.394073   0.017060   0.008826
75%     0.499345   0.534064   0.197912
max     0.919071   1.181992   1.127869
A          B          C
count  20.000000  20.000000  20.000000
mean    0.576878   0.612413   0.508070
std     0.442297   0.466569   0.575018
min    -0.223874  -0.175731  -0.570708
25%     0.261960   0.352292  -0.031961
50%     0.572278   0.619588   0.677623
75%     0.926991   0.916843   0.836673
max     1.418941   1.608790   1.320450


Voici quelques représentations graphiques des données. D'abord variable par variable :

In [4]:
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')

Out[4]:
<matplotlib.text.Text at 0x10edc6850>
In [5]:
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')

Out[5]:
[<matplotlib.lines.Line2D at 0x10f106b10>]

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 :

In [6]:
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')

Out[6]:
[<matplotlib.lines.Line2D at 0x10f3ca510>]

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 :

In [7]:
fig = plt.figure()
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')

Out[7]:
<matplotlib.text.Text at 0x10f3eded0>

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 :

In [8]:
df = pd.concat([jeu1,jeu2],ignore_index=True)
print 'df :', df.shape, '\n', df.head()

df : (40, 3)
A         B         C
0  0.448471 -0.571576  0.192196
1  0.504769 -0.047861 -1.713761
2 -0.531144  0.511938  0.024027
3  0.780623 -0.305116  0.507643
4 -0.703031  0.600445 -1.058459

In [9]:
df.describe()

Out[9]:
A B C
count 40.000000 40.000000 40.000000
mean 0.390655 0.337737 0.213275
std 0.491479 0.562774 0.646706
min -0.703031 -0.626166 -1.713761
25% 0.116848 -0.135448 -0.105172
50% 0.436283 0.453357 0.148563
75% 0.733606 0.664721 0.740718
max 1.418941 1.608790 1.320450

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é.

In [10]:
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'

U : dimensions : (40, 40)
[[-0.08016633  0.24954351  0.02017695 ..., -0.05241184  0.05128136
0.03639094]
[-0.41447148 -0.06072098 -0.28138753 ...,  0.07287957  0.19822456
0.10318354]
[-0.0906032  -0.13580931  0.260539   ..., -0.23240311 -0.29179327
-0.15403108]
...,
[ 0.14304098  0.10155904  0.15281051 ...,  0.94529937 -0.06059889
-0.03173766]
[ 0.2183952  -0.01492947  0.26114769 ..., -0.06450452  0.90758327
-0.04911905]
[ 0.10973732 -0.01300494  0.14325562 ..., -0.03291673 -0.04838817
0.97419441]]
lambda :
[ 4.30736035  3.3720568   2.85642057]
V :dimensions : (3, 3)
[[ 0.31667818  0.37972899  0.86920701]
[ 0.2611152  -0.91586366  0.3049797 ]
[-0.91188474 -0.13038275  0.38918704]]



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 :

In [11]:
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 :

In [12]:
fig2 = plt.figure(figsize=(7,7))
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")
plt.title('Eigenvectors')
plt.xlabel('A')
plt.ylabel('B')

0.316678175078 0.379728987521 0.86920701186
0.261115200802 -0.915863655166 0.304979699417
-0.911884743451 -0.130382748793 0.389187041765

Out[12]:
<matplotlib.text.Text at 0x10f99d150>

Calculons les composantes principales :

In [13]:
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']

Out[13]:
CP1 CP2 CP3
0 -0.345305 0.841475 0.057634
1 -1.785278 -0.204755 -0.803761
2 -0.390261 -0.457957 0.744209
3 0.135251 0.780368 -0.157225
4 -1.351989 -0.914035 0.468121

Représentons la première composante principale et comparez aux graphiques obtenus plus haut :

In [14]:
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')

Out[14]:
<matplotlib.text.Text at 0x10faaff90>

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 :

In [15]:
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')

Out[15]:
<matplotlib.text.Text at 0x10fb8ed10>

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.

In [16]:
n = df_centre.shape[0]
inertie_tot = ((df_centre**2).sum()).sum()/n
print 'inertie totale : ', inertie_tot

inertie totale :  0.952081466702


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 :

In [17]:
sum(lamb**2)/n

Out[17]:
0.95208146670208627

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$ :

In [18]:
lamb[0]**2/n

Out[18]:
0.46383382895239167

et pour la proportion de la variance expliquée on a :

In [19]:
lamb[0]**2/sum(lamb**2)

Out[19]:
0.48717871860174428

ç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 :

In [20]:
df.var()*(n-1)/sum(lamb**2)

Out[20]:
A    0.247366
B    0.324338
C    0.428296
dtype: float64

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$ :

In [21]:
sum(lamb[0:2]**2)/n

Out[21]:
0.74810300547304964

Et pour la part de la variance expliquée par les deux premières composantes principales :

In [22]:
sum(lamb[0:2]**2)/sum(lamb**2)

Out[22]:
0.78575524431160559

Le module sklearn contient une fonction pour faire l'ACP qui est plus commode à utiliser :

In [23]:
from sklearn.decomposition import PCA

In [24]:
sklearn_pca = PCA()
df_transf = sklearn_pca.fit_transform(df)  # pas nécessaire de centrer les colonnes ; renvoie les comp. princ.
df_transf

Out[24]:
array([[ -3.45305277e-01,   8.41474884e-01,   5.76338533e-02],
[ -1.78527803e+00,  -2.04754598e-01,  -8.03761127e-01],
[ -3.90260652e-01,  -4.57956709e-01,   7.44208966e-01],
[  1.35251013e-01,   7.80368183e-01,  -1.57224703e-01],
[ -1.35198866e+00,  -9.14035467e-01,   4.68120908e-01],
[ -1.58123088e-01,  -4.93031204e-01,  -5.53732237e-02],
[ -7.79662074e-01,   5.25884660e-01,   3.57905320e-01],
[ -1.04420679e+00,   3.74937007e-01,  -3.72479624e-01],
[ -3.19475478e-02,   2.50459291e-01,  -5.18852961e-01],
[ -3.45311645e-01,   1.74335913e-01,  -2.68721175e-01],
[  1.70569210e-03,   7.52858419e-02,  -5.82177865e-02],
[ -3.45913370e-01,   8.42914004e-01,   4.64442037e-01],
[ -1.01629811e+00,   1.57253031e-01,   4.75694505e-01],
[  1.03268398e-02,  -2.68019172e-01,  -4.32064607e-01],
[ -3.14677504e-01,   5.87058902e-01,   5.53877464e-02],
[ -2.39181819e-01,   6.53838162e-01,   6.74896004e-02],
[ -7.57342821e-02,  -1.00957411e+00,   3.95901630e-01],
[ -3.16997962e-01,  -6.12431675e-01,   8.40849227e-01],
[  6.61122622e-01,   3.34512519e-01,   6.11331915e-01],
[ -6.57788958e-01,   6.22153141e-01,  -5.43478502e-02],
[  1.83397356e-01,  -7.56656353e-01,   3.81665005e-02],
[  2.55522816e-01,  -9.20453253e-01,  -7.37688720e-01],
[  5.23776081e-01,  -6.93475270e-02,   2.28581490e-01],
[  4.26286410e-01,   6.90541490e-01,  -8.29219006e-02],
[  9.74369911e-01,  -1.84278943e-01,  -5.38437268e-01],
[  8.73938910e-01,  -9.35985928e-01,  -6.28592183e-01],
[  2.57938913e-02,  -5.77337505e-01,  -2.34635551e-01],
[  4.75641513e-01,  -5.54669965e-01,   1.18733274e-01],
[ -1.39889441e-01,   3.92792760e-01,  -6.27526058e-01],
[ -6.82282854e-01,  -6.27055826e-01,   2.57412479e-01],
[  3.14638093e-01,  -1.37945964e-01,   4.47601030e-01],
[  1.07046519e+00,   3.46991039e-01,   1.82159392e-01],
[  1.27378806e+00,   2.57973169e-01,  -2.65552766e-01],
[  6.17086738e-01,   4.15113164e-01,   2.23771333e-01],
[ -3.30250573e-03,  -3.21071025e-01,  -4.11775006e-01],
[  9.19137516e-01,   1.35074212e-01,  -7.69774364e-01],
[ -7.47612157e-01,   3.37377391e-01,  -6.09081331e-01],
[  6.16129049e-01,   3.42462850e-01,   4.36491072e-01],
[  9.40706840e-01,  -5.03430060e-02,   7.45947639e-01],
[  4.72678185e-01,  -4.38533890e-02,   4.09198289e-01]])
In [25]:
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')

Out[25]:
<matplotlib.text.Text at 0x111a9bf90>

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 :

In [26]:
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)

Inertie d'une, deux, trois composantes principales :
[ 0.46383383  0.74810301  0.95208147]
valeurs singulieres : [ 4.30736035  3.3720568   2.85642057]


# 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

1. Importez les données.
2. Familiarisez-vous avec les données. Vérifiez qu'il y a bien le même nombre d'observations pour chaque espèce.
3. 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.
4. Effectuez une ACP et représentez graphiquement les résultats.
5. Calculez l'inertie totale et l'inertie des données projetées.
In [ ]: