Travaux pratiques - Analyse en composantes principales

(correspond à 1 séance de TP)

_images/jupyter_logo.png

Cahier Jupyter

Références externes utiles :

L’objectif de cette séance de TP est de présenter l’utilisation des fonctionnalités de Scikit-learn concernant l’ACP, ainsi que de contribuer à une meilleure compréhension de cette méthode d’analyse de données. Après une brève présentation de la classe correspondante de Scikit-learn, nous examinons d’abord des données générées de façon contrôlée afin de mieux comprendre l’impact de la distribution des données sur les résultats de l’ACP. Deux exemples du cours, basés sur des données réelles, sont ensuite repris et examinés.

L’ACP dans Scikit-learn

Dans Scikit-learn, l’analyse en composantes principales (ACP) est mise en œuvre dans la classe sklearn.decomposition.PCA (voir http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html). Les variantes et quelques exemples sont présentés sur http://scikit-learn.org/stable/modules/decomposition.html.

Pour définir l’analyse on appelle PCA(n_components=None, copy=True, whiten=False, svd_solver='auto', tol=0.0, iterated_power='auto', random_state=None). Les paramètres d’appel sont (pour des explications plus détaillées voir le lien ci-dessus) :

  • n_components : nombre de composantes à conserver (par défaut toutes, c’est à dire min(n_samples, n_features) où n_samples est le nombre de lignes de la matrice de données et n_features le nombre de colonnes).

  • copy : si False, la matrice de données est écrasée par les données transformées (par défaut True).

  • whiten : si True, les données transformées sont modifiées pour que les projections sur les axes principaux présentent une variance unitaire (par défaut False).

  • svd_solver : solveur à employer parmi 'auto', 'full', 'arpack', 'randomized' (par défaut 'auto').

  • tol, iterated_power, random_state : paramètres optionnels pour solveurs spécifiques.

Les attributs accessibles sont les suivants :

  • explained_variance_ : array [n_components] dans lequel on trouve les valeurs propres en ordre décroissant.

  • components_ : array [n_components, n_features] dans lequel chaque ligne correspond à un vecteur propre ; l’ordre est celui des valeurs propres.

  • explained_variance_ratio_ : array [n_components] avec les valeurs propres exprimées en pourcentage de la variance expliquée.

  • mean_ : array [n_features] contenant les moyennes des variables calculées lors de l’ACP.

  • n_components_ : nombre estimé de composantes.

  • noise_variance_ : variance du bruit estimée suivant l’approche de l’ACP probabiliste (voir lien ci-dessus).

Les méthodes qui peuvent être employées :

  • fit(X[, y]) : trouver les axes principaux et les valeurs propres associés à partir de la matrice de données X.

  • fit_transform(X[, y]) : trouver les axes principaux et les valeurs propres associés, appliquer la projection sur les axes principaux (avec stockage des résultats dans la matrice de données X si copy=False et dans une autre matrice si copy=True).

  • get_covariance() : calculer les covariances à partir de la matrice de données.

  • get_params([deep]) : retourner les paramètres de l’estimateur.

  • get_precision() : calculer la matrice de précision des données.

  • inverse_transform(X[, y]) : re-projeter les données transformées dans l’espace initial.

  • score(X[, y]) : valeur moyenne de la log-vraisemblance pour toutes les observations.

  • score_samples(X) : log-vraisemblance pour chaque observation.

  • set_params(**params) : donner des valeurs aux paramètres de l’estimateur.

  • transform(X[, y]) : appliquer la projection sur les axes principaux (avec stockage des résultats dans la matrice de données X si copy=False et dans une autre matrice si copy=True).

ACP de données générées

Démarrez par la génération d’un ensemble de 500 vecteurs dans l’espace 3D, suivant une loi normale (de moyenne nulle et de variance unitaire), ensuite visualisez les points générés :

import numpy as np
# Génération de données selon une loi normale tridimensionnelle
rndn3d = np.random.randn(500,3)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Affichage du nuage de points
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(rndn3d[:,0], rndn3d[:,1], rndn3d[:,2])
plt.title("Données initiales")
plt.show()

Appliquez ensuite l’ACP à ces données :

from sklearn.decomposition import PCA

pca = PCA(n_components=3)
pca.fit(rndn3d)

print("Pourcentage de variance expliquée : ")
print(pca.explained_variance_ratio_)
print("Composantes principales : ")
print(pca.components_)

Chaque ligne de pca.components_ représente le vecteur unitaire qui donne la direction d’un axe factoriel. Le vecteur pca.components_[i,:] correspond à la valeur propre pca.explained_variance_ratio_[i].

Question :

Générez un autre tableau de données en utilisant de nouveau np.random.randn(500,3) et appliquez l’ACP à ces nouvelles données. Comparez les résultats avec ceux obtenus sur rndn3d. Que constatez-vous ? Expliquez.

Correction :

On constate d’abord que les trois valeurs propres sont chaque fois très proches entre elles, ce qui indique que le nuage d’observations a une forme relativement sphérique. Ensuite, en examinant les vecteurs propres correspondants, on constate que d’un ensemble de données générées à un autre les directions des vecteurs propres changent, donc ces vecteurs n’ont pas de pertinence particulière pour décrire les données. Ce qui est cohérent avec les faibles différences entre valeurs propres successives.

Appliquez maintenant une déformation et une rotation dans l’espace tridimensionnel au nuage des observations de rndn3d, visualisez le résultat :

s1 = np.array([[3,0,0],[0,1,0],[0,0,0.2]])  # matrice de déformation
r1 = np.array([[0.36,0.48,-0.8],[-0.8,0.6,0],[0.48,0.64,0.6]])  # matrice de rotation
rndef = rndn3d.dot(s1).dot(r1)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(rndef[:,0], rndef[:,1], rndef[:,2])
plt.title("Données déformées")
plt.show()

Appliquez une ACP à ces observations :

pca = PCA(n_components=3)
pca.fit(rndef)

print("Pourcentage de variance expliquée : ")
print(pca.explained_variance_ratio_)
print("Composantes principales : ")
print(pca.components_)

Question :

Que constatez-vous par rapport à l’application directe de l’ACP sur les données de rndn3d ? Expliquez. Peut-on lier les rapports entre les valeurs propres obtenues aux rapports entre les éléments de la diagonale de la matrice de déformation s1 ?

Correction :

La décroissance des valeurs propres est ici très rapide, le premier axe factoriel explique 90% de la variance. Les rapports entre valeurs propres correspondent aux carrés des rapports entre les éléments de la diagonale de la matrice de déformation. En effet, les valeurs propres sont des variances alors que les éléments de la diagonale de la matrice de déformation ont ici une signification d’écart-types.

Question :

Générez un autre ensemble de données en utilisant les mêmes matrices de transformation et appliquez l’ACP à ces nouvelles données. Comparez les résultats avec ceux obtenus sur rndef. Que constatez-vous ? Expliquez.

Correction :

Les valeurs propres obtenues sur le second ensemble de données sont très proches de celles obtenues sur rndef. De plus, entre les deux ensembles de données, les vecteurs propres correspondant aux valeurs propres de même ordre sont très similaires. Les valeurs et vecteurs propres caractérisent les (fortes) déformations appliquées aux nuages d’observations générés et sont donc relativement stables par rapport à un changement d’échantillon.

ACP des données sur le sommeil des mammifères

L’ACP sera d’abord appliquée aux données concernant le sommeil des mammifères, issues de [AC76] et déjà vues dans le cours. Pour cela, il est d’abord nécessaire d’enregistrer ces données dans le répertoire de travail. Dans une fenêtre terminal entrez (ou télécharger le fichier via ce lien):

wget -nc http://cedric.cnam.fr/~crucianm/src/mammals.csv

Vous pouvez examiner le contenu de ce fichier format .csv. Les données doivent être ensuite lues dans Python :

mammals = np.loadtxt('mammals.csv', delimiter=';', usecols=[1,2,3,4,5,6,7,8,9,10], skiprows=1)
print(mammals[:2,:])
noms = np.genfromtxt('mammals.csv', dtype='str', delimiter=';', usecols=[0], skip_header=1)
print(noms)

Question :

Appliquez une ACP aux données de mammals et affichez les valeurs propres. Comparez les deux premières aux valeurs vues en cours. Que constatez-vous ? Pourquoi ?

Correction :

pca = PCA()
pca.fit(mammals)

print(pca.explained_variance_ratio_)
# array([  9.61492649e-01,   3.35026815e-02,   4.90760735e-03,
#          8.11547630e-05,   1.21589728e-05,   2.12907093e-06,
#          8.08341160e-07,   4.40101019e-07,   3.18913491e-07,
#          5.23943276e-08])

La première valeur propre obtenue explique la quasi-totalité de la variance, alors que dans le cours elle expliquait moins de la moitié. L’ACP directe n’est pas une ACP normée, les variances des variables « poids » (du corps, du cerveau) sont comparativement très élevées par rapport aux variances des autres variables et dominent très largement dans la définition de la direction du premier axe factoriel. Pour avoir une analyse utile il faut centrer et réduire les variables.

Question :

Appliquez une transformation de centrage et réduction aux variables décrivant les données de mammals (regardez le centrage et la réduction des variables (« standardization ») avec scikit-learn) et conservez le résultat dans mammalsCR. Appliquez de nouveau l’ACP, de préférence en utilisant une nouvelle instance pcaCR pour conserver pca. Comment ont évolué les valeurs propres les plus grandes ? Affichez le graphique de décroissance des valeurs propres.

Correction :

from sklearn import preprocessing

mammalsCR = preprocessing.scale(mammals)
pcaCR = PCA()
pcaCR.fit(mammalsCR)

ratios = pcaCR.explained_variance_ratio_
print(ratios)
# array([ 0.48017524,  0.22078825,  0.12605711,  0.06488067,  0.04739057,
#         0.02607844,  0.01797744,  0.00920976,  0.00501433,  0.00242818])
plt.bar(range(len(ratios)), ratios)
plt.xticks(range(len(ratios)))
plt.xlabel("Composante principale")
plt.ylabel("% de variance expliquée")
plt.show()

La plus grande valeur propre a maintenant une valeur qui correspond au pourcentage vu en cours.

Nous pouvons maintenant obtenir et afficher les projections des observations sur les deux premiers axes factoriels. Rappelons que mammalsCR contient les données transformées (variables centrées et réduites) et que pcaCR est l’instance de PCA utilisée pour analyser mammalsCR ; mammalsCR et pcaCR ont été obtenues en réponse à la question ci-dessus.

mNt = pcaCR.transform(mammalsCR)
print(mNt[:2,:])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111)
for i in range(len(noms)):
        x,y = mNt[i,0], mNt[i,1]
        ax.scatter(x,y)
        ax.text(x,y,noms[i])

plt.show()

Question :

Réalisez l’affichage des projections sur les 3 premiers axes principaux ; la possibilité de rotation du graphique facilite la lecture des étiquettes textuelles.

Correction :

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
for i in range(len(noms)):
        x,y,z = mNt[i,0],mNt[i,1],mNt[i,2]
        ax.scatter(x,y,z)
        ax.text(x,y,z,noms[i])

plt.show()

Question :

Les résultats de l’analyse du nuage des variables ne sont pas directement disponibles. Comment pouvez-vous les obtenir ? Affichez la projection du nuage des variables sur les deux premiers axes factoriels.

Correction :

Il suffit pour cela d’utiliser les relations de transition entre les deux analyses (voir le cours).

mammalsCRTranspose = mammalsCR.transpose()
irlambdas = 1/np.sqrt(pcaCR.explained_variance_)
mirlambdas = np.diagflat(irlambdas)
projectionsVars = (mammalsCRTranspose.dot(mNt)).dot(mirlambdas)/62
nomsVars = np.genfromtxt('mammals.csv', dtype='str', delimiter=';', skip_header=0, max_rows=1)

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
cercle = plt.Circle((0, 0), 1, color='b', fill=False)
ax.add_artist(cercle)
for i in range(len(nomsVars)-1):
        x,y = projectionsVars[i,0], projectionsVars[i,1]
        ax.scatter(x,y)
        ax.arrow(0, 0, x, y, shape='full', lw=1, length_includes_head=True)
        ax.text(x, y, nomsVars[i+1])

plt.show()

ACP sur les données « feuilles »

Nous appliquerons maintenant l’ACP sur les données concernant des feuilles d’arbres, issues de https://archive.ics.uci.edu/ml/datasets/Leaf. Pour cela, il est d’abord nécessaire d’enregistrer ces données dans le répertoire de travail. Dans une fenêtre terminal entrez (ou téléchargez depuis cette adresse:

wget -nc http://cedric.cnam.fr/~crucianm/src/leaf.csv

Examinez le fichier téléchargé et regardez sur le site indiqué les explications concernant le format des données.

Appliquez l’ACP aux données et affichez la décroissance des valeurs propres :

leaf = np.loadtxt('leaf.csv', delimiter=',')
print(leaf[:2,:])

from sklearn import preprocessing
leafN = preprocessing.scale(leaf[:,2:])
pca = PCA()
pca.fit(leafN)

ratios = pca.explained_variance_ratio_

plt.bar(range(len(ratios)), pca.explained_variance_ratio_)
plt.xticks(range(len(ratios)))
plt.xlabel("Composante principale")
plt.ylabel("% de variance expliquée")
plt.show()

Question :

Affichez les projections des données sur les 3 premiers axes principaux.

Correction :

leafNt = pca.transform(leafN)
leafNt[:2,:]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(leafNt[:,0],leafNt[:,1],leafNt[:,2])
plt.show()

Question :

Affichez les projections des données sur les 3 premiers axes principaux en utilisant l’étiquette de classe pour donner une couleur aux marques. L’étiquette de classe se trouve (sous forme numérique) dans la première colonne de leaf. Il est nécessaire d’ajouter dans l’appel ax.scatter un argument pour indiquer que la couleur est obtenue à partir de données, c = ....

Correction :

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(leafNt[:,0],leafNt[:,1],leafNt[:,2],c=leaf[:,0])
plt.show()
[AC76]

Allison T., D. Cicchetti. Sleep in mammals: ecological and constitutional correlates. Science, 194: 732–734 (1976).