python - principales - pca como funciona
Análisis de componentes principales(PCA) en Python (9)
Tengo una matriz (26424 x 144) y quiero realizar una PCA con Python. Sin embargo, no hay un lugar particular en la web que explique cómo lograr esta tarea (hay algunos sitios que solo hacen PCA de acuerdo con los suyos, no hay una forma generalizada de hacerlo que pueda encontrar). Cualquiera con algún tipo de ayuda lo hará genial.
Además de todas las otras respuestas, aquí hay un código para trazar el biplot usando sklearn y matplotlib .
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.decomposition import PCA
import pandas as pd
from sklearn.preprocessing import StandardScaler
iris = datasets.load_iris()
X = iris.data
y = iris.target
#In general a good idea is to scale the data
scaler = StandardScaler()
scaler.fit(X)
X=scaler.transform(X)
pca = PCA()
x_new = pca.fit_transform(X)
def myplot(score,coeff,labels=None):
xs = score[:,0]
ys = score[:,1]
n = coeff.shape[0]
scalex = 1.0/(xs.max() - xs.min())
scaley = 1.0/(ys.max() - ys.min())
plt.scatter(xs * scalex,ys * scaley, c = y)
for i in range(n):
plt.arrow(0, 0, coeff[i,0], coeff[i,1],color = ''r'',alpha = 0.5)
if labels is None:
plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, "Var"+str(i+1), color = ''g'', ha = ''center'', va = ''center'')
else:
plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, labels[i], color = ''g'', ha = ''center'', va = ''center'')
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.xlabel("PC{}".format(1))
plt.ylabel("PC{}".format(2))
plt.grid()
#Call the function. Use only the 2 PCs.
myplot(x_new[:,0:2],np.transpose(pca.components_[0:2, :]))
plt.show()
Aquí hay opciones de scikit-learn. Con ambos métodos, se utilizó StandardScaler porque PCA se efectúa por escala
Método 1: haga que scikit-learn elija la cantidad mínima de componentes principales, de modo que al menos x% (90% en el ejemplo a continuación) de la varianza se conserve.
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
iris = load_iris()
# mean-centers and auto-scales the data
standardizedData = StandardScaler().fit_transform(iris.data)
pca = PCA(.90)
principalComponents = pca.fit_transform(X = standardizedData)
# To get how many principal components was chosen
print(pca.n_components_)
Método 2: elija la cantidad de componentes principales (en este caso, se eligió 2)
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
iris = load_iris()
standardizedData = StandardScaler().fit_transform(iris.data)
pca = PCA(n_components=2)
principalComponents = pca.fit_transform(X = standardizedData)
# to get how much variance was retained
print(pca.explained_variance_ratio_.sum())
Fuente: https://towardsdatascience.com/pca-using-python-scikit-learn-e653f8989e60
Este es un trabajo para numpy .
Y aquí hay un tutorial que demuestra cómo el análisis de componentes principales se puede hacer usando los numpy incorporados de numpy como mean,cov,double,cumsum,dot,linalg,array,rank .
http://glowingpython.blogspot.sg/2011/07/principal-component-analysis-with-numpy.html
Observe que scipy también tiene una explicación larga aquí - https://github.com/scikit-learn/scikit-learn/blob/babe4a5d0637ca172d47e1dfdd2f6f3c3ecb28db/scikits/learn/utils/extmath.py#L105
con la biblioteca scikit-learn con más ejemplos de código - https://github.com/scikit-learn/scikit-learn/blob/babe4a5d0637ca172d47e1dfdd2f6f3c3ecb28db/scikits/learn/utils/extmath.py#L105
He hecho un pequeño guión para comparar las diferentes PCA que aparecieron como respuesta aquí:
import numpy as np
from scipy.linalg import svd
shape = (26424, 144)
repeat = 20
pca_components = 2
data = np.array(np.random.randint(255, size=shape)).astype(''float64'')
# data normalization
# data.dot(data.T)
# (U, s, Va) = svd(data, full_matrices=False)
# data = data / s[0]
from fbpca import diffsnorm
from timeit import default_timer as timer
from scipy.linalg import svd
start = timer()
for i in range(repeat):
(U, s, Va) = svd(data, full_matrices=False)
time = timer() - start
err = diffsnorm(data, U, s, Va)
print(''svd time: %.3fms, accuracy: %E'' % (time*1000/repeat, err))
from matplotlib.mlab import PCA
start = timer()
_pca = PCA(data)
for i in range(repeat):
U = _pca.project(data)
time = timer() - start
err = diffsnorm(data, U, _pca.fracs, _pca.Wt)
print(''matplotlib PCA time: %.3fms, accuracy: %E'' % (time*1000/repeat, err))
from fbpca import pca
start = timer()
for i in range(repeat):
(U, s, Va) = pca(data, pca_components, True)
time = timer() - start
err = diffsnorm(data, U, s, Va)
print(''facebook pca time: %.3fms, accuracy: %E'' % (time*1000/repeat, err))
from sklearn.decomposition import PCA
start = timer()
_pca = PCA(n_components = pca_components)
_pca.fit(data)
for i in range(repeat):
U = _pca.transform(data)
time = timer() - start
err = diffsnorm(data, U, _pca.explained_variance_, _pca.components_)
print(''sklearn PCA time: %.3fms, accuracy: %E'' % (time*1000/repeat, err))
start = timer()
for i in range(repeat):
(U, s, Va) = pca_mark(data, pca_components)
time = timer() - start
err = diffsnorm(data, U, s, Va.T)
print(''pca by Mark time: %.3fms, accuracy: %E'' % (time*1000/repeat, err))
start = timer()
for i in range(repeat):
(U, s, Va) = pca_doug(data, pca_components)
time = timer() - start
err = diffsnorm(data, U, s[:pca_components], Va.T)
print(''pca by doug time: %.3fms, accuracy: %E'' % (time*1000/repeat, err))
pca_mark es el pca en la respuesta de Marcos .
pca_doug es la respuesta de pca en doug .
Aquí hay un ejemplo de salida (pero el resultado depende en gran medida del tamaño de los datos y de los componentes pca, por lo que recomendaría ejecutar su propia prueba con sus propios datos. Además, el pca de Facebook está optimizado para datos normalizados, por lo que será más rápido y más preciso en ese caso):
svd time: 3212.228ms, accuracy: 1.907320E-10
matplotlib PCA time: 879.210ms, accuracy: 2.478853E+05
facebook pca time: 485.483ms, accuracy: 1.260335E+04
sklearn PCA time: 169.832ms, accuracy: 7.469847E+07
pca by Mark time: 293.758ms, accuracy: 1.713129E+02
pca by doug time: 300.326ms, accuracy: 1.707492E+02
Otra PCA Python usando numpy. La misma idea que @doug pero esa no se ejecutó.
from numpy import array, dot, mean, std, empty, argsort
from numpy.linalg import eigh, solve
from numpy.random import randn
from matplotlib.pyplot import subplots, show
def cov(data):
"""
Covariance matrix
note: specifically for mean-centered data
note: numpy''s `cov` uses N-1 as normalization
"""
return dot(X.T, X) / X.shape[0]
# N = data.shape[1]
# C = empty((N, N))
# for j in range(N):
# C[j, j] = mean(data[:, j] * data[:, j])
# for k in range(j + 1, N):
# C[j, k] = C[k, j] = mean(data[:, j] * data[:, k])
# return C
def pca(data, pc_count = None):
"""
Principal component analysis using eigenvalues
note: this mean-centers and auto-scales the data (in-place)
"""
data -= mean(data, 0)
data /= std(data, 0)
C = cov(data)
E, V = eigh(C)
key = argsort(E)[::-1][:pc_count]
E, V = E[key], V[:, key]
U = dot(data, V) # used to be dot(V.T, data.T).T
return U, E, V
""" test data """
data = array([randn(8) for k in range(150)])
data[:50, 2:4] += 5
data[50:, 2:5] += 5
""" visualize """
trans = pca(data, 3)[0]
fig, (ax1, ax2) = subplots(1, 2)
ax1.scatter(data[:50, 0], data[:50, 1], c = ''r'')
ax1.scatter(data[50:, 0], data[50:, 1], c = ''b'')
ax2.scatter(trans[:50, 0], trans[:50, 1], c = ''r'')
ax2.scatter(trans[50:, 0], trans[50:, 1], c = ''b'')
show()
Lo cual produce lo mismo que el mucho más corto
from sklearn.decomposition import PCA
def pca2(data, pc_count = None):
return PCA(n_components = 4).fit_transform(data)
Tal como lo entiendo, usar valores propios (primera vía) es mejor para datos de alta dimensión y menos muestras, mientras que el uso de la descomposición de valores singulares es mejor si tiene más muestras que dimensiones.
Por el bien def plot_pca(data): funcionará, es necesario reemplazar las líneas
data_resc, data_orig = PCA(data)
ax1.plot(data_resc[:, 0], data_resc[:, 1], ''.'', mfc=clr1, mec=clr1)
con líneas
newData, data_resc, data_orig = PCA(data)
ax1.plot(newData[:, 0], newData[:, 1], ''.'', mfc=clr1, mec=clr1)
Publiqué mi respuesta a pesar de que ya se ha aceptado otra respuesta; la respuesta aceptada se basa en una matplotlib.org/api/… ; adicionalmente, esta función en desuso se basa en la Descomposición de Valor Singular (SVD), que (aunque es perfectamente válida) es mucho más intensiva en memoria y procesador de las dos técnicas generales para calcular PCA. Esto es particularmente relevante aquí debido al tamaño de la matriz de datos en el PO. Usando una PCA basada en covarianza, la matriz utilizada en el flujo de cálculo es solo 144 x 144 , en lugar de 26424 x 144 (las dimensiones de la matriz de datos original).
Aquí hay una implementación simple de trabajo de PCA usando el módulo linalg de SciPy . Debido a que esta implementación primero calcula la matriz de covarianza y luego realiza todos los cálculos posteriores en esta matriz, usa mucha menos memoria que la PCA basada en SVD.
(El módulo linalg en NumPy también se puede usar sin cambios en el código a continuación, aparte de la declaración de importación, que sería desde numpy import linalg como LA ).
Los dos pasos clave en esta implementación de PCA son:
calcular la matriz de covarianza ; y
tomando los eivenvectors y valores propios de esta matriz de cov
En la función siguiente, el parámetro dims_rescaled_data se refiere a la cantidad deseada de dimensiones en la matriz de datos redimensionados ; este parámetro tiene un valor predeterminado de solo dos dimensiones, pero el código siguiente no está limitado a dos, pero podría ser cualquier valor menor que el número de columna de la matriz de datos original.
def PCA(data, dims_rescaled_data=2):
"""
returns: data transformed in 2 dims/columns + regenerated original data
pass in: data as 2D NumPy array
"""
import numpy as NP
from scipy import linalg as LA
m, n = data.shape
# mean center the data
data -= data.mean(axis=0)
# calculate the covariance matrix
R = NP.cov(data, rowvar=False)
# calculate eigenvectors & eigenvalues of the covariance matrix
# use ''eigh'' rather than ''eig'' since R is symmetric,
# the performance gain is substantial
evals, evecs = LA.eigh(R)
# sort eigenvalue in decreasing order
idx = NP.argsort(evals)[::-1]
evecs = evecs[:,idx]
# sort eigenvectors according to same index
evals = evals[idx]
# select the first n eigenvectors (n is desired dimension
# of rescaled data array, or dims_rescaled_data)
evecs = evecs[:, :dims_rescaled_data]
# carry out the transformation on the data using eigenvectors
# and return the re-scaled data, eigenvalues, and eigenvectors
return NP.dot(evecs.T, data.T).T, evals, evecs
def test_PCA(data, dims_rescaled_data=2):
''''''
test by attempting to recover original data array from
the eigenvectors of its covariance matrix & comparing that
''recovered'' array with the original data
''''''
_ , _ , eigenvectors = PCA(data, dim_rescaled_data=2)
data_recovered = NP.dot(eigenvectors, m).T
data_recovered += data_recovered.mean(axis=0)
assert NP.allclose(data, data_recovered)
def plot_pca(data):
from matplotlib import pyplot as MPL
clr1 = ''#2026B2''
fig = MPL.figure()
ax1 = fig.add_subplot(111)
data_resc, data_orig = PCA(data)
ax1.plot(data_resc[:, 0], data_resc[:, 1], ''.'', mfc=clr1, mec=clr1)
MPL.show()
>>> # iris, probably the most widely used reference data set in ML
>>> df = "~/iris.csv"
>>> data = NP.loadtxt(df, delimiter='','')
>>> # remove class labels
>>> data = data[:,:-1]
>>> plot_pca(data)
El siguiente diagrama es una representación visual de esta función PCA en los datos del iris. Como puede ver, una transformación en 2D separa limpiamente la clase I de la clase II y la clase III (pero no la clase II de la clase III, que de hecho requiere otra dimensión).
Puede encontrar una función PCA en el módulo matplotlib:
import numpy as np
from matplotlib.mlab import PCA
data = np.array(np.random.randint(10,size=(10,3)))
results = PCA(data)
los resultados almacenarán los diversos parámetros de la PCA. Es de la parte mlab de matplotlib, que es la capa de compatibilidad con la sintaxis de MATLAB
EDIT: en el blog nextgenetics , encontré una maravillosa demostración de cómo realizar y mostrar una PCA con el módulo matplotlib mlab, ¡diviértanse y revisen ese blog!
ACTUALIZACIÓN: matplotlib.mlab.PCA desde la versión 2.2 (2018-03-06) de hecho en deprecated .
La biblioteca matplotlib.mlab.PCA (utilizada en esta respuesta ) no está en desuso. Entonces, para toda la gente que llega aquí a través de Google, publicaré un ejemplo de trabajo completo probado con Python 2.7.
Utilice el siguiente código con cuidado ya que utiliza una biblioteca ahora en desuso.
from matplotlib.mlab import PCA
import numpy
data = numpy.array( [[3,2,5], [-2,1,6], [-1,0,4], [4,3,4], [10,-5,-6]] )
pca = PCA(data)
Ahora en `pca.Y ''está la matriz de datos original en términos de los principales vectores de base de componentes. Más detalles sobre el objeto PCA se pueden encontrar here .
>>> pca.Y
array([[ 0.67629162, -0.49384752, 0.14489202],
[ 1.26314784, 0.60164795, 0.02858026],
[ 0.64937611, 0.69057287, -0.06833576],
[ 0.60697227, -0.90088738, -0.11194732],
[-3.19578784, 0.10251408, 0.00681079]])
Puede usar matplotlib.pyplot para dibujar estos datos, solo para convencerse de que PCA produce "buenos" resultados. La lista de names solo se usa para anotar nuestros cinco vectores.
import matplotlib.pyplot
names = [ "A", "B", "C", "D", "E" ]
matplotlib.pyplot.scatter(pca.Y[:,0], pca.Y[:,1])
for label, x, y in zip(names, pca.Y[:,0], pca.Y[:,1]):
matplotlib.pyplot.annotate( label, xy=(x, y), xytext=(-2, 2), textcoords=''offset points'', ha=''right'', va=''bottom'' )
matplotlib.pyplot.show()
Al mirar nuestros vectores originales, veremos que los datos [0] ("A") y los datos [3] ("D") son bastante similares, como lo son los datos [1] ("B") y los datos [2] (" DO"). Esto se refleja en el gráfico 2D de nuestros datos transformados PCA.
