Saltar al contenido

Análisis de componentes principales (PCA) en Python

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).

ingrese la descripción de la imagen aquí

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.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *