Travaux pratiques - Analyse en composantes principales

(correspond à 1 séance de TP)

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

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 ?

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.

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 ?

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.

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.

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.

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.

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


[AC76]

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