In [36]:
import numpy as np
import pandas as pa
import matplotlib.pyplot as plt
from sklearn import decomposition
from sklearn import preprocessing
In [37]:
data = pa.read_csv('https://www.labri.fr/~zemmari/datasets/eleves.csv', decimal='.', index_col=0, header=0)
data.head()
Out[37]:
Math Sciences Francais Latin Musique
Prenom
Jean 6.0 6.0 5.0 5.5 8
Aline 8.0 8.0 8.0 8.0 9
Annie 6.0 7.0 11.0 9.5 11
Monique 14.5 14.5 15.5 15.0 8
Didier 14.0 14.0 12.0 12.0 10
In [56]:
pca_cols = ['Math', 'Sciences', 'Francais', 'Latin', 'Musique']
data_pca = data[pca_cols]
X = data_pca.values
X
n = X.shape[0]
p = X.shape[1]
In [101]:
corr = np.corrcoef(X, rowvar=0)

corr
Out[101]:
array([[ 1.        ,  0.98253573,  0.22673193,  0.49049826,  0.01118384],
       [ 0.98253573,  1.        ,  0.39669324,  0.6339855 ,  0.00630993],
       [ 0.22673193,  0.39669324,  1.        ,  0.95613107,  0.03803599],
       [ 0.49049826,  0.6339855 ,  0.95613107,  1.        ,  0.08860159],
       [ 0.01118384,  0.00630993,  0.03803599,  0.08860159,  1.        ]])
In [57]:
std_scaler = preprocessing.StandardScaler()
Z = std_scaler.fit_transform(X)
Z
Out[57]:
array([[-1.08650256, -1.28173989, -1.50366322, -1.62520902, -1.01904933],
       [-0.4938648 , -0.61300603, -0.63985669, -0.72231512, -0.67936622],
       [-1.08650256, -0.94737296,  0.22394984, -0.18057878,  0.        ],
       [ 1.43220792,  1.560379  ,  1.51965963,  1.8057878 , -1.01904933],
       [ 1.28404848,  1.39319553,  0.51188535,  0.72231512, -0.33968311],
       [ 0.39509184,  0.05572782, -1.35969546, -1.08347268,  0.67936622],
       [-1.234662  , -0.94737296,  1.08775637,  0.54173634, -0.33968311],
       [ 0.9877296 ,  0.89164514, -0.49588893, -0.18057878,  0.33968311],
       [-0.19754592, -0.11145564,  0.65585311,  0.72231512,  2.37778177]])
In [58]:
#On peut vérifier que nos données sont bien centrées réduites :
print(np.mean(Z,axis=0))
print(np.std(Z,axis=0))
[  1.81953218e-16  -2.20502629e-16   2.34380416e-16  -1.23358114e-17
   4.93432455e-17]
[ 1.  1.  1.  1.  1.]
In [70]:
nb_composantes = 5
pca = decomposition.PCA(n_components=nb_composantes)
coord_fact = pca.fit_transform(Z)
In [71]:
valeurs_propres = pca.explained_variance_
valeurs_propres
Out[71]:
array([  3.21954472e+00,   1.29451627e+00,   1.10603333e+00,
         4.42926059e-03,   4.76420198e-04])
In [72]:
var_cumulee = valeurs_propres.cumsum()
var_cumulee
Out[72]:
array([ 3.21954472,  4.51406099,  5.62009432,  5.62452358,  5.625     ])
In [73]:
#On peut également afficher le pourcentage de la variance expliquées par chaque valeur propre:
variance_par_vp = pca.explained_variance_ratio_*100
variance_par_vp
Out[73]:
array([  5.72363506e+01,   2.30136226e+01,   1.96628148e+01,
         7.87424105e-02,   8.46969241e-03])
In [74]:
#Et la variance accumulée : 
variances_acc = variance_par_vp.cumsum()
variances_acc
Out[74]:
array([  57.23635057,   80.24997313,   99.9127879 ,   99.99153031,  100.        ])
In [75]:
plt.plot(np.arange(1,p+1),valeurs_propres) 
plt.title("Scree plot") 
plt.ylabel("Eigen values") 
plt.xlabel("Factor number") 
plt.show() 
In [76]:
#cumul de variance expliquée 
plt.plot(np.arange(1,p+1),np.cumsum(pca.explained_variance_ratio_)) 
plt.title("Explained variance vs. # of factors") 
plt.ylabel("Cumsum explained variance ratio") 
plt.xlabel("Factor number") 
plt.show() 
In [77]:
#Les coordonnées factorielles des individus est données par l'instriction 
#coord_fact = pca.fit_transform(Z) --> Voir la cellule correspondante
coord_fact
Out[77]:
array([[-2.78571626, -0.6764554 , -0.73680069,  0.04824298, -0.0332307 ],
       [-1.26249031, -0.33033839, -0.55495038,  0.03576028,  0.01091186],
       [-1.01674674,  1.01982862, -0.28808844, -0.06237732,  0.00743856],
       [ 3.12214541, -0.16584307, -1.14426837,  0.10012136,  0.00532702],
       [ 1.95510007, -0.78790499, -0.18921203, -0.09873943, -0.01278293],
       [-0.9477522 , -1.20140468,  1.14011052,  0.01857612,  0.04027278],
       [-0.32496763,  1.75473914, -0.90951345, -0.04309011,  0.01406045],
       [ 0.63734131, -1.12985744,  0.69185868, -0.05383744, -0.01383015],
       [ 0.62308636,  1.51723621,  1.99086415,  0.05534356, -0.01816687]])
In [82]:
#contributions aux axes 
ctr = coord_fact**2 
for j in range(p):     
    ctr[:,j] = ctr[:,j]/(n*valeurs_propres[j])      
print(pa.DataFrame({'id':data.index,'CTR_1':ctr[:,0],'CTR_2':ctr[:,1]}))
      CTR_1     CTR_2        id
0  0.267816  0.039276      Jean
1  0.055007  0.009366     Aline
2  0.035677  0.089270     Annie
3  0.336410  0.002361   Monique
4  0.131917  0.053284    Didier
5  0.030999  0.123888     Andre
6  0.003645  0.264287    Pierre
7  0.014019  0.109571  Brigitte
8  0.013399  0.197586   Evelyne
In [84]:
#On ne garde que les individus ayant une contribution > à la contribution moyenne 1/n = 1/9
1/9
#Donc pour CP1, il s'git de Jean, Monique et Didier 
# et pour CP2, il s'agit de Andre, Pierre, Brigitte (~) et Evelyne
#==> On attends l'interprétation pour les variables pour réinterpréter ces résultats
Out[84]:
0.1111111111111111
In [85]:
#On récupère d'abords les vecteurs propres :
vect_propres = pca.components_
print(vect_propres)
[[ 0.47639106  0.53024454  0.44812425  0.53807342  0.03941992]
 [-0.53263677 -0.40159625  0.56963554  0.37055187  0.305328  ]
 [ 0.15475344  0.09364365 -0.22763607 -0.10934894  0.95052978]
 [ 0.30298298 -0.51679947 -0.47749743  0.64155887 -0.03896199]
 [ 0.61123017 -0.53076393  0.44144992 -0.38679883  0.01399912]]
In [87]:
#racine carrée des valeurs propres 
sqrt_eigval = np.sqrt(valeurs_propres)
In [90]:
#corrélation des variables avec les axes 
corvar = np.zeros((p,p)) 
for k in range(p):     
        corvar[:,k] = pca.components_[k,:] * sqrt_eigval[k]      
#afficher la matrice des corrélations variables x facteurs     
print(corvar) 
#Les variables sont maintenant en ligne, les facteurs en colonne : 
[[  8.54792763e-01  -6.06017128e-01   1.62751286e-01   2.01643322e-02
    1.33413526e-02]
 [  9.51422546e-01  -4.56923404e-01   9.84832657e-02  -3.43943954e-02
   -1.15850118e-02]
 [  8.04073365e-01   6.48113152e-01  -2.39400574e-01  -3.17787388e-02
    9.63555027e-03]
 [  9.65469972e-01   4.21602099e-01  -1.15000219e-01   4.26974688e-02
   -8.44267810e-03]
 [  7.07315146e-02   3.47392466e-01   9.99654295e-01  -2.59302509e-03
    3.05559567e-04]]
In [91]:
#Si l’on s’en tient spécifiquement aux deux premiers facteurs : 

#on affiche pour les deux premiers axes 
print(pa.DataFrame({'id':data.columns,'COR_1':corvar[:,0],'COR_2':corvar[:,1]})) 
      COR_1     COR_2        id
0  0.854793 -0.606017      Math
1  0.951423 -0.456923  Sciences
2  0.804073  0.648113  Francais
3  0.965470  0.421602     Latin
4  0.070732  0.347392   Musique
In [95]:
#Contribution des variables aux axes (CTR). 
#contributions  
ctrvar = corvar**2
for k in range(p):     
    ctrvar[:,k] = ctrvar[:,k]/valeurs_propres[k] 
 #on n'affiche que pour les deux premiers axes 
print(pa.DataFrame({'id':data.columns,'CTR_1':ctrvar[:,0],'CTR_2':ctrvar[:,1]})) 
      CTR_1     CTR_2        id
0  0.226948  0.283702      Math
1  0.281159  0.161280  Sciences
2  0.200815  0.324485  Francais
3  0.289523  0.137309     Latin
4  0.001554  0.093225   Musique
In [96]:
#On ne garde alors que les variables ayant une ctr > à la contribution moyenne 
1/p
Out[96]:
0.2
In [ ]:
#Donc pour CP1: on garde Math, Sciences, Français et latin
#et pour CP2 : Math, Français