PCA en numpy et sklearn produit des résultats différents

21

Suis-je mal comprendre quelque chose. C'est mon code

en utilisant sklearn

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import decomposition
from sklearn import datasets
from sklearn.preprocessing import StandardScaler

pca = decomposition.PCA(n_components=3)

x = np.array([
        [0.387,4878, 5.42],
        [0.723,12104,5.25],
        [1,12756,5.52],
        [1.524,6787,3.94],
    ])
pca.fit_transform(x)

Production:

array([[ -4.25324997e+03,  -8.41288672e-01,  -8.37858943e-03],
   [  2.97275001e+03,  -1.25977271e-01,   1.82476780e-01],
   [  3.62475003e+03,  -1.56843494e-01,  -1.65224286e-01],
   [ -2.34425007e+03,   1.12410944e+00,  -8.87390454e-03]])

Utilisation de méthodes numpy

x_std = StandardScaler().fit_transform(x)
cov = np.cov(x_std.T)
ev , eig = np.linalg.eig(cov)
a = eig.dot(x_std.T)

Production

array([[ 0.06406894,  0.94063993, -1.62373172],
   [-0.35357757,  0.7509653 ,  0.63365168],
   [ 0.29312477,  0.6710958 ,  1.11766206],
   [-0.00361615, -2.36270102, -0.12758202]])
I have kept all 3 components but it doesnt seem to allow me to retain my original data.

Puis-je savoir pourquoi en est-il ainsi?

Si je veux récupérer ma matrice d'origine, que dois-je faire?

aceminer
la source
Votre code numpy est faux à mon humble avis (également, il utilise Xce qui n'est pas défini). Revérifiez vos mathématiques .
Anony-Mousse -Reinstate Monica
J'utilise le bloc-notes ipython, donc je ne pouvais copier que par cellules. Mes maths sont faux? Quelle partie @ Anony-Mousse
aceminer
@ Anony-Mousse Oui j'ai réalisé mon erreur mais elle ne correspond toujours pas
aceminer
@aceminer Je suis curieux de savoir pourquoi vous calculez la matrice de covariance de x_std.T, pas x_std?
Evgeni Nabokov
@EvgeniNabokov ça fait trop longtemps. Sry je ne me souviens pas déjà
aceminer

Réponses:

21

La différence est parce decomposition.PCAque ne standardise pas vos variables avant de faire PCA, alors que dans votre calcul manuel vous appelez StandardScalerpour faire la standardisation. Vous observez donc cette différence: l' ACP sur la corrélation ou la covariance?

Si vous remplacez

pca.fit_transform(x)

avec

x_std = StandardScaler().fit_transform(x)
pca.fit_transform(x_std)

vous obtiendrez le même résultat qu'avec le calcul manuel ...

... mais uniquement jusqu'à l'ordre des PC. En effet, lorsque vous exécutez

ev , eig = np.linalg.eig(cov)

vous obtenez des valeurs propres pas nécessairement dans l'ordre décroissant. Je reçois

array([ 0.07168571,  2.49382602,  1.43448827])

Vous voudrez donc les commander manuellement. Sklearn fait cela pour vous.


En ce qui concerne la reconstruction des variables d'origine, voir Comment inverser l'ACP et reconstruire les variables d'origine à partir de plusieurs composants principaux?

amibe dit réintégrer Monica
la source
Je voudrais juste vérifier. Est-il vraiment nécessaire de standardiser la matrice par son écart type? J'ai vu des exemples où ils ne le font pas
aceminer
Ce n'est pas nécessaire , c'est juste une façon de le faire. Voir le lien que j'ai mis dans le premier paragraphe: stats.stackexchange.com/questions/53 - il s'agit vraiment de cette question. Si vous standardisez, vous effectuez l'ACP sur les corrélations. Si vous ne le faites pas, vous faites l'APC sur les covariances.
Amoeba dit Reinstate Monica
9

Voici une belle implémentation avec discussion et explication de PCA en python. Cette implémentation conduit au même résultat que le scikit PCA. Ceci est un autre indicateur que votre PCA est erroné.

import numpy as np
from scipy import linalg as LA

x = np.array([
        [0.387,4878, 5.42],
        [0.723,12104,5.25],
        [1,12756,5.52],
        [1.524,6787,3.94],
    ])

#centering the data
x -= np.mean(x, axis = 0)  

cov = np.cov(x, rowvar = False)

evals , evecs = LA.eigh(cov)

vous devez trier les valeurs propres (et les vecteurs propres en conséquence) en ordre décroissant

idx = np.argsort(evals)[::-1]
evecs = evecs[:,idx]
evals = evals[idx]

a = np.dot(x, evecs) 

Généralement, je vous recommande de vérifier votre code en implémentant un exemple simple (aussi simple que possible) et en calculant à la main les bons résultats (et résultats intermédiaires). Cela vous aide à identifier le problème.

Nikolas Rieble
la source
1
J'adore cette réponse. Cela a résolu mon problème!
Jinhua Wang