Solución:
Publiqué mi respuesta a pesar de que ya se ha aceptado otra respuesta; la respuesta aceptada se basa en una función obsoleta; Además, esta función obsoleta se basa en Valor singular de descomposición (SVD), que (aunque perfectamente válida) es la 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 OP. Usando PCA basado en covarianza, la matriz utilizada en el flujo de cálculo es simplemente 144 x 144, en vez de 26424 x 144 (las dimensiones de la matriz de datos original).
Aquí hay una implementación funcional simple de PCA usando el linalg módulo de Ciencia. Debido a que esta implementación calcula primero la matriz de covarianza y luego realiza todos los cálculos posteriores en esta matriz, utiliza mucha menos memoria que el PCA basado en SVD.
(el módulo de linalg en NumPy también se puede utilizar sin cambios en el código siguiente, aparte de la declaración de importación, que sería de numpy import linalg como LA.)
Los dos pasos clave en esta implementación de PCA son:
-
calculando el Matriz de covarianza; y
-
tomando el eivenvectores Y valores propios de esta cov matriz
En la función siguiente, el parámetro dims_rescaled_data se refiere al número deseado de dimensiones en el reescalado matriz de datos; este parámetro tiene un valor predeterminado de solo dos dimensiones, pero el código siguiente no se limita a dos, pero podría ser alguna 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 gráfico siguiente es una representación visual de esta función de PCA en los datos del iris. Como puede ver, una transformación 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 del PCA. Es de la parte mlab de matplotlib, que es la capa de compatibilidad con la sintaxis de MATLAB.
EDITAR: en el blog nextgenetics encontré una maravillosa demostración de cómo realizar y mostrar un PCA con el módulo matplotlib mlab, ¡diviértete y revisa ese blog!
Otro Python PCA usando numpy. La misma idea que @doug pero esa no funcionó.
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(X):
"""
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 que 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)
Según tengo entendido, usar valores propios (primera forma) es mejor para datos de alta dimensión y menos muestras, mientras que usar la descomposición de valores singulares es mejor si tiene más muestras que dimensiones.