Résumé: Ce calepin introduit l'utilisation de la librairie scikit-learn pour l'exploration statistique. Ceci est illustré par des exemples de mise en oeuvre de l'(ACP sur des données "jouet" puis sur des images élémentaires de caractères et enfin sur des données économiques sous une la forme particulière d'un cube ou tableauà trois indices.
Scikit-learn vs. RL'objectif de ce tutoriel est d'introduire l'utilisation de la librairie scikit-learn de Python pour l'exploration de données multidimensionnelles. Seule l'utilisation directe des fonctions d'exploration est abordée. Se pose rapidement une question: quand utiliser scikit-learn de Python plutôt que R plus complet et plus simple d'emploi?
Le choix repose sur les points suivants:
array de numpy chargés en mémoire et donc de taille limitée par la RAM de l'ordinateur; de façon analogue R charge en RAM des objets de type data.frame.scikit-learn (0.18) ne reconnaît pas (ou pas encore ?) la classe DataFrame de pandas; scikit-learn utilise la classe array de numpy. C'est un problème pour la gestion de variables qualitatives complexes. Une variable binaire est simplement remplacée par un codage $(0,1)$ mais, en présence de plusieurs modalités, traiter celles-ci comme des entiers n'a pas de sens statistique et remplacer une variable qualitative par l'ensemble des indicatrices (dummy variables $(0,1)$) de ses modalités complique l'interprétation statistique. scikit-learn sont aussi efficaces (e.g. $k$-means), voire beaucoup plus efficaces pour des données volumineuses car utilisent implicitement les capacités de parallélisation.scikit-learn se fait sous contrôle d'un groupe qui en contrôle la pertinence et l'efficacité. En conséquences:
scikit-learn pour mettre au point une chaîne de traitements (pipe line) opérationnelle de l'extraction à une analyse privilégiant la prévision brute à l'interprétation et pour des données quantitatives ou rendues quantitatives ("vectorisation" de corpus de textes).scikit-learnLa communauté qui développe cette librairie est très active, elle évolue rapidement. Ne pas hésiter à consulter la documentation pour des compléments. Voici une sélection de ses principales fonctionnalités.
scikit-learnL'objectif est d'illustrer la mise en oeuvre de l'analyse en composantes principales. Consulter la documentation et ses nombreux exemples pour plus de détails sur l'utilisation de scikit-learn.
La librairie scikit-learn a principalement été conçu en vue de l'analyse de signaux. Aussi, de nombreuses options de l'ACP ne sont pas disponibles, notamment les graphiques usuels (biplot, cercle des corrélations...). En revanche des résultats sont liés à la version probabiliste de l'ACP sous hypothèse d'une distribution gaussienne multidimensionnelle des données.
Attention, l'ACP est évidemment centrée mais par réduite. L'option n'est pas prévue et les variables doivent être réduites (fonction sklearn.preprocessing.scale) avant si c'est nécessaire. L'attribut transform désigne les composantes principales, sous-entendu: transformation par réduction de la dimension; n_components fixe le nombre de composantes retenues, par défaut toutes; l'attribut components_ contient les n_components vecteurs propres mais non normalisés, c'est-à-dire de norme carrée la valeur propre associée, et donc à utiliser pour représenter les variables.
D'autres versions d'analyse en composantes principales sont proposées dans Scikit-learn: kernel PCA, sparse PCA, ICA...
Plusieurs jeux de données élémentaires sont utilisés donyt celui "jouet" déjà vu en R afin de bien comprendre les sorties proposées par la fonction disponible. L'autre ensemble de données est un problème classique et simplifié de reconnaissance de caractères qui est inclus dans la librairie scikit-learn.
Les données sont celles de l'exemple introduction à l'ACP: les notes en maths, français, physique et anglais de 9 lycéens virtuels. L'objectif est de contrôler les résultats en les comparant avec ceux obtenus avec R.
C'est une façon générique de procéder à l'approche d'un nouveau logiciel ou de fonctionnalités inconnues: traiter des donées triviales dont les résultats de l'analyse sont parfaitement maîtrisés.
In [ ]:
# Construire la matrice de notes
import pandas as pd
note=[[6,6,5,5.5],[8,8,8,8],[6,7,11,9.5],[14.5,14.5,15.5,15],
[14,14,12,12.5],[11,10,5.5,7],[5.5,7,14,11.5],[13,12.5,8.5,9.5],
[9,9.5,12.5,12]]
dat=pd.DataFrame(note,index=["jean","alai","anni","moni","didi","andr","pier","brig","evel"],
columns=["Math","Phys","Fran","Angl"])
dat
In [ ]:
# Importation des fonctions
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
import numpy as np
Attention Les valeurs singulières sont celles de la décomposition de la matrice centrée par rapport aux métriques usuelles: $(\bar{X}, I_p, I_n)$ alors que le diviseur de la variance est celui d'une estimation sans biais: $(n-1)$.
Contrairement à beaucoup de logiciels, l'ACP de scikit-learn n'est pas réduite.
In [ ]:
pca = PCA()
pca.fit(dat).explained_variance_
In [ ]:
pca.singular_values_
Les valeurs singulières associées à l'ACP sont celles de $(\bar{X}, I_p, \frac{1}{n-1}I_n)$
In [ ]:
pca.singular_values_/np.sqrt(8)
Pour retrouver les valeurs propres de l'ACP à partir des valeurs singulières de la matrice centrée:
In [ ]:
(pca.singular_values_/np.sqrt(8))**2
In [ ]:
pca.components_.T
In [ ]:
pca.transform(dat)
Q Comparer avec les résultats obtenus en R.
Tous les autres résultats (contributions, cossinus carrés, corrélations variables facteurs...) et surtout les graphes (éboulis, plans factoriels...) sont à construire car aucune fonction n'est disponible comme dans FactoMineR. C'est partièlement fait dans le jeu de données suivant et complété (biplot) dans les calepins plus completes des cas d'usage.
Il s'agit d'explorer les données issues de la pixellisation de tracés de caractères dont les procédés d'obtention et prétraitement sont décrits sur le site de l'UCI (Lichman, 2013). Les chiffres ont été saisies sur des tablettes à l'intérieur de cadres de résolution $500\times 500$. Des procédures de normalisation, ré-échantillonnage spatial puis de lissage ont été appliquées. Chaque caractère apparaît finalement discrétisé sous la forme d'une matrice $8\times 8$ de pixels à 16 niveaux de gris et identifié par un label. Les données sont archivées sous la forme d'une matrice ou tableau à trois indices. Elles sont également archivées après vectorisation des images sous la forme d'une matrice à $p=64$ colonnes.
L'étude du même type de données, mais nettement plus complexes (MNIST): 60 000 caractères représentés par des images de 784 pixels (26 $\times$ 26) fait l'objet d'un autre calepin.
In [ ]:
# Importations
import matplotlib.pyplot as plt
from sklearn import datasets
%matplotlib inline
# les données présentes dnas la librairie
digits = datasets.load_digits()
# Contenu et mode d'obtention
print(digits)
In [ ]:
# Dimensions
digits.images.shape
In [ ]:
# Sous forme d'un cube d'images 1797 x 8x8
print(digits.images)
In [ ]:
# Sous forme d'une matrice 1797 x 64
print(digits.data)
In [ ]:
# Label réel de chaque caractère
print(digits.target)
Voici un aperçu des empilements des images à décrire puis ensuite en principe à discriminer:
In [ ]:
images_and_labels = list(zip(digits.images,
digits.target))
for index, (image, label) in enumerate(images_and_labels[:8]):
plt.subplot(2, 4, index + 1)
plt.axis('off')
plt.imshow(image, cmap=plt.cm.gray_r, interpolation='nearest')
plt.title('Chiffres: %i' % label)
In [ ]:
from sklearn.decomposition import PCA
X=digits.data
y=digits.target
target_name=[0,1,2,3,4,5,6,7,8,9]
# définition de la commande
pca = PCA()
# Estimation, calcul des composantes principales
C = pca.fit(X).transform(X)
# Décroissance de la variance expliquée
plt.plot(pca.explained_variance_ratio_)
plt.show()
Diagramme boîte des premières composantes principales.
In [ ]:
plt.boxplot(C[:,0:20])
plt.show()
Q Quelle dimension retenir en principe?
Représentation des caractères dans le premier plan principal.
La représentation des variables (pixels) et le biplot n'ont pas grand intérêt pour ces données.
In [ ]:
plt.scatter(C[:,0], C[:,1], c=y, label=target_name)
plt.show()
Le même graphique avec une légende mais moins de couleurs.
In [ ]:
# attention aux indentations
plt.figure()
for c, i, target_name in zip("rgbcmykrgb",[0,1,2,3,4,5,6,7,8,9], target_name):
plt.scatter(C[y == i,0], C[y == i,1], c=c, label=target_name)
plt.legend()
plt.title("ACP Digits")
plt.show()
Graphique en trois dimensions.
In [ ]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(1, figsize=(8, 6))
ax = Axes3D(fig, elev=-150, azim=110)
ax.scatter(C[:, 0], C[:, 1], C[:, 2], c=y, cmap=plt.cm.Paired)
ax.set_title("ACP: trois premieres composantes")
ax.set_xlabel("Comp1")
ax.w_xaxis.set_ticklabels([])
ax.set_ylabel("Comp2")
ax.w_yaxis.set_ticklabels([])
ax.set_zlabel("Comp3")
ax.w_zaxis.set_ticklabels([])
plt.show()
L'objectif de cette section est l'exploration de données socio-économiques plus complexes. La principale spécificité de ces données est de se présenter sous la forme d'un cube de données ou tableau à trois entrées: le numéro de ligne, le numéro de variable et l'année d'observation de cette variable. Après une description classique, la mise en oeuvre de l'analyse en composantes principales avec python nécessite un effort particulier afin de produire les graphes adaptés à la structure particulière des données.
Les données sont issues de l'Observatoire de l'OCDE. Pour chaque pays membre et pour chacune des années 1975, 1977, 1979, 1981, on connaît les valeurs prises par les variables suivantes qui sont toutes des \emph{taux}~:
Elles sont disponibles dans le fichier: ocdeR.dat.
Les mêmes variables sont donc observées, sur les mêmes pays ou individus à quatre dates différentes. Plusieurs stratégies d'analyse sont possibles (tableau moyen, tableaux concaténés, meilleur compromis ou double ACP). La plus adaptée pour ces données est de considérer les observations des variables pour chacun des individus: pays $\times$ années.
In [ ]:
# Importaiton des principals librairies et
# Affichage des graphiques dans le notebook
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
In [ ]:
ocde=pd.read_table("Data/ocdeR.dat",sep='\s+',index_col=0)
ocde.head()
In [ ]:
ocde.mean()
In [ ]:
ocde["CNRJ"].hist(bins=20)
plt.show()
In [ ]:
from pandas.plotting import scatter_matrix
scatter_matrix(ocde, alpha=0.2, figsize=(15, 15), diagonal='kde')
plt.show()
Q Quel est le graphique ci-dessous? Que représentent les blocs dagonaux? Que dire des structures de corrélation?
Chaque pays étant observé 4 fois, la principale difficulté technique est de faire apparaître cette structure chronologique dans les graphique afin d'illustrer la dynamique économique de la période considérée.
Q Justifier la nécessité de réduire.
Q Pourqoi toutes les variables sont des taux?
In [ ]:
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
# réduction
ocdeS=scale(ocde)
pca = PCA()
cpOcde = pca.fit_transform(ocdeS)
# Eboulis
plt.plot(pca.explained_variance_ratio_)
plt.show()
In [ ]:
plt.boxplot(cpOcde)
plt.show()
In [ ]:
coord1=pca.components_[0]*np.sqrt(pca.explained_variance_[0])
coord2=pca.components_[1]*np.sqrt(pca.explained_variance_[1])
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(1, 1, 1)
for i, j, nom in zip(coord1,coord2, ocde.columns):
plt.text(i, j, nom)
plt.arrow(0,0,i,j,color='black')
plt.axis((-1.2,1.2,-1.2,1.2))
# cercle
c=plt.Circle((0,0), radius=1, color='gray', fill=False)
ax.add_patch(c)
plt.show()
Q Interpréter chacun des deux premiers axes.
Exo représenter le plan (2,3) et interpréter le 3ème axe.
In [ ]:
plt.figure(figsize=(10,6))
for i, j, nom in zip(cpOcde[:,0], cpOcde[:,1], ocde.index):
# color = int(i/4)
plt.text(i, j, nom ,color="blue")
plt.axis((-5,7,-4,4))
plt.show()
In [ ]:
import matplotlib.patheffects as PathEffects
comp_0 = 0
comp_1 = 1
cmap = plt.get_cmap("tab20")
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(1,1,1)
for i,k in enumerate(np.arange(0,cpOcde.shape[0],4)):
country =ocde.index[k]
xs = cpOcde[k:k+4,comp_0]
ys = cpOcde[k:k+4, comp_1]
ax.plot(xs,ys, color=cmap(i), marker=".", markersize=15)
txt = ax.text(xs[-4], ys[-4], country, horizontalalignment="left", verticalalignment="top",
color=cmap(i), fontweight="bold", fontsize=15)
# Add black line around text
#txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='black')])
ax.set_xlabel("PC%d" %comp_0, fontsize=20)
ax.set_ylabel("PC%d" %comp_1, fontsize=20)
plt.tight_layout()
plt.show()
Q Analyser les évolutions des économies des différents pays. Les remplacer dans la période considérée.