Travaux pratiques - Analyse factorielle discriminante

(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’AFD, 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 l’application de l’AFD à des données générées de façon contrôlée, sans ou avec une réduction de dimension préalable par ACP. Deux exemples du cours, basés sur des données réelles, sont ensuite repris et analysés.

L’AFD dans Scikit-learn

Dans Scikit-learn, l’analyse factorielle discriminante (AFD) est mise en œuvre dans la classe LinearDiscriminantAnalysis). Le guide utilisateur fournit quelques explications supplémentaires.

Pour définir l’analyse, on appelle LinearDiscriminantAnalysis(solver='svd', shrinkage=None, priors=None, n_components=None, store_covariance=False, tol=0.0001).

Note

Attention, la classe sklearn.lda.LDA est obsolète (deprecated). Il ne faut plus l’utiliser.

Les paramètres d’appel sont (pour des explications plus détaillées voir le lien ci-dessus) :

  • solver : solveur employé pour obtenir les axes ; 'svd' traite directement la matrice de données, sans calculer au préalable la matrice des covariances empiriques ; 'lsrq' est la solution des moindres carrés ; 'eigen' calcule d’abord la matrice des covariances empiriques (par défaut, solver = 'svd').

  • shrinkage : si absent, la matrice des covariances empiriques est directement utilisée ; si 'auto', choix automatique du paramètre suivant le lemme Ledoit-Wolf [LWH04].

  • priors : optionnel, array n_classes avec les a priori des différentes classes.

  • n_components : nombre de composantes à conserver (par défaut toutes, c’est à dire < n_classes - 1)).

  • store_covariance : optionnel, si True demande le calcul de la maytrice des covariances empiriques.

  • tol : optionnel, seuil à utiliser pour l’estimation du rang dans le solveur SVD.

Les attributs accessibles sont les suivants :

  • coef_ : array, unidimensionnel n_features ou bidimensionnel n_classes x n_features, poids des variables, éventuellement par classe.

  • intercept_ : array n_features, intercept.

  • covariance_ : array n_features x n_features, matrice des covariances empiriques (partagée par toutes les classes).

  • explained_variance_ratio_ : array n_components, pourcentage de variance expliquée par chaque composante (ou facteur discriminant).

  • means_ : array n_classes x n_features, centres de gravité des classes.

  • priors_ : array n_classes avec les a priori des différentes classes (leur somme est égale à 1).

  • scalings_ : array rank x n_classes - 1, transformation des variables dans l’espace des centres de gravité des classes.

  • xbar_ : array n_features, centre de gravité global.

  • classes_ : array n_classes, étiquettes des classes.

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

  • decision_function(X) : scores de confiance pour les données.

  • fit(X, y[, store_covariance, tol]) : construction du modèle à partir des données X et des étiquettes de classe y.

  • fit_transform(X[, y]) : construction du modèle et ensuite projection des données sur l’espace discriminant.

  • get_params([deep]) : lire les valeurs des paramètres de l’estimateur employé.

  • predict(X) : prédire les étiquettes de classe pour les données de la matrice X (par ex. données de test).

  • predict_log_proba(X) : estimer les log-probabilités.

  • predict_proba(X) : estimer les probabilités.

  • score(X, y[, sample_weight]) : taux de bon classement des données X, avec les étiquettes désirées dans y.

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

  • transform(X) : projeter les données sur l’espace discriminant.

AFD de données générées

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

import numpy as np    # si pas encore fait
s1 = np.array([[3,0,0],[0,1,0],[0,0,0.01]])  # 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

# On génère deux nuages de points déformés et tournés
rndn3d1 = np.random.randn(500,3)
rndef1 = rndn3d1.dot(s1).dot(r1)
rndn3d2 = np.random.randn(500,3)
# Le deuxième nuage de point est translaté selon l'axe z
rndef2 = rndn3d2.dot(s1).dot(r1) + [0, 0, 1]
rndef = np.concatenate((rndef1, rndef2))
print(rndef.shape)

Question :

Expliquez l’addition qui permet de créer rndef2 : à quoi sert cette addition, quel est le résultat de l’addition entre deux array de dimension différente ?

Correction :

Cette addition sert à décaler un des deux nuages (représenté par les données de rndef2) sur l’axe z et ainsi à séparer les données des deux classes. Pour l’addition, le vecteur [0, 0, 1] est (implicitement) utilisé pour produire une matrice de 500 lignes (comme rndn3d2.dot(s1).dot(r1)) dont chaque ligne est égale à ce vecteur.

Maintenant, générez des étiquettes de classe correspondantes : 1 pour les données de rndef1, 2 pour celles de rndef2 et ensuite visualisez les données :

lcls1 = np.ones(500)
lcls2 = 2 * np.ones(500)
lcls = np.concatenate((lcls1, lcls2))
print(lcls)

Cela indique que les 500 premiers points sont de classe 1 et les 500 points suivants sont de classe 2.

Affichons maintenant le nuage de points avec matplotlib :

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

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

Appliquez d’abord une analyse discriminante à ces données, projetez-les sur l’espace discriminant (de quelle dimension est-il ?) et ensuite visualisez les résultats :

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

lda = LinearDiscriminantAnalysis()
lda.fit(rndef,lcls)

L’AFD a ainsi été entraînée sur les données. Nous pouvons maintenant obtenir les données projetées grâce à la méthode transform :

rndt = lda.transform(rndef)
print(rndt.shape)

Ce qui produit un vecteur de dimension (1000, 1), c’est-à-dire que l’AFD a projeté nos observations sur un unique axe discriminant.

Affichons désormais le résultat de cette projection :

plt.plot(rndt, lcls, 'r+')
plt.xlabel("Axe discriminant")
plt.ylabel("Classe")
plt.yticks([1, 2])
plt.show()

Essayez de réduire la dimension des données avec une ACP en conservant les projections sur les 2 premiers axes principaux, appliquez ensuite l’AFD et visualisez les résultats :

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(rndef)

Projection des observations puis affichage du graphe :

rndp = pca.transform(rndef)
print(rndp.shape)
# (1000, 2)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(rndp[:,0], rndp[:,1], c=lcls)
plt.show()

On applique ensuite une AFD par-dessus le résultat de cette projection :

lda = LinearDiscriminantAnalysis()
lda.fit(rndp,lcls)
rndpt = lda.transform(rndp)

plt.plot(rndpt, lcls, 'r+')
plt.xlabel("Axe discriminant")
plt.ylabel("Classe")
plt.yticks([1, 2])
plt.show()

Question :

Expliquez ce que vous constatez.

Correction :

L’ACP préalable a supprimé la direction qui permettait de séparer les deux classes car dans cette direction la variance du nuage global était la plus faible. En conséquence, les classes ne sont plus séparables pour les projections des observations sur ce premier plan factoriel.

ACP et AFD des données « textures »

Ces données correspondent à 5500 observations qui sont des pixels extraits d’images de textures et décrits par 40 variables (descripteurs visuels). Chaque observation appartient à une des 11 classes de textures ; chaque classe est représentée par 500 pixels. Les données sont issues de https://www.elen.ucl.ac.be/neural-nets/Research/Projects/ELENA/databases/REAL/texture/

wget -nc http://cedric.cnam.fr/~crucianm/src/texture.dat

Examinez le fichier téléchargé et regardez sur le site indiqué ci-dessus les explications concernant le format des données. Le jeu de données comporte 5500 observations comportant 41 colonnes : 40 attributs visuels décrivant le type de texture et 1 attribut concernant à la catégorie de la texture.

Question :

Lisez ces données dans un array en utilisant np.loadtxt. Appliquez une ACP et visualisez les projections sur les 3 premiers axes principaux. Appliquez également une AFD (sur les données initiales, non sur les projections issues de l’ACP) et visualisez les projections sur les 3 premiers axes discriminants.

Correction :

textures = np.loadtxt('texture.dat')

pca = PCA()
pca.fit(textures[:,:40])
texturesp = pca.transform(textures[:,:40])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(texturesp[:,0], texturesp[:,1], texturesp[:,2], c=textures[:,40])
plt.title("Projection après ACP")
plt.show()

lda.fit(textures[:,:40],textures[:,40])
texturest = lda.transform(textures[:,:40])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(texturest[:,0], texturest[:,1], texturest[:,2], c=textures[:,40])
plt.title("Projection après AFD")
plt.show()

Question :

Faites l’extraction d’une partie des observations avant d’appliquer l’AFD. Réalisez ensuite la prédiction (étape décisionnelle de l’AFD) pour ces observations et affichez l’erreur de prédiction. Pour l’extraction des observations vous pouvez employer (entre autres) :

from sklearn.model_selection import train_test_split
train, test = train_test_split(textures, test_size = 0.2)

Correction :

from sklearn.model_selection import train_test_split
train, test = train_test_split(textures, test_size = 0.2)
print(train.shape)
# (4400, 41)
print(test.shape)
# (1100, 41)

lda.fit(train[:,:40],train[:,40])
lda.score(test[:,:40],test[:,40])
# 0.99818181818181817

AFD des données « feuilles »

Nous appliquerons enfin l’AFD sur les données concernant des feuilles d’arbres, issues de https://archive.ics.uci.edu/ml/datasets/Leaf. Si le données ne se trouvent pas déjà dans le répertoire de travail, entrez dans une fenêtre terminal :

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

Examinez une nouvelle fois le fichier téléchargé et regardez sur le site indiqué ci-dessus les explications concernant le format des données.

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

lda = LinearDiscriminantAnalysis()
lda.fit(leaf[:,2:],leaf[:,0])

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

Question :

Comparez les résultats avec ceux de l’ACP.

Correction :

Pour comparer il est nécessaire de reprendre l’affichage des projections des observations sur les trois axes principaux en utilisant pour chaque point la couleur correspondant à la classe, comme nous l’avons fait ci-dessus pour les projections obtenues par AFD. En comparant les deux graphiques on constate que l’AFD sépare un peu mieux les classes, mais celles-ci restent difficiles à déceler car leur nombre est très important (40 !) et donc les nuances de couleur correspondant à des classes différentes sont parfois très proches. L’utilisation de la couleur pour distinguer les classes a atteint ses limites.

Question :

Faites l’extraction d’une partie des observations avant d’appliquer l’AFD. Réalisez ensuite la prédiction (étape décisionnelle de l’AFD) pour ces observations et affichez l’erreur de prédiction.

Correction :

from sklearn.model_selection import train_test_split
train, test = train_test_split(leaf, test_size = 0.2)
lda.fit(train[:,2:],train[:,0])
lda.score(test[:,2:],test[:,0])
# 0.8970588235294118

[LWH04]

Ledoit O., Wolf M. Honey, I Shrunk the Sample Covariance Matrix. Journal of Portfolio Management 30(4), 110-119, 2004.