Principal component analysis (PCA) is a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components. The number of principal components is less than or equal to the number of original variables. This transformation is defined in such a way that the first principal component has the largest possible variance (that is, accounts for as much of the variability in the data as possible), and each succeeding component in turn has the highest variance possible under the constraint that it is orthogonal to the preceding components. The resulting vectors are an uncorrelated orthogonal basis set. The principal components are orthogonal because they are the eigenvectors of the covariance matrix, which is symmetric. PCA is sensitive to the relative scaling of the original variables.
In [ ]:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
In [ ]:
def gaussiana2D(prom, sigma, cov, n):
np.random.seed(142)
return np.random.multivariate_normal(np.array([prom, prom]), np.array([[sigma, cov], [cov, sigma]]), n)
In [ ]:
datosNoCorrel = gaussiana2D(prom=50, sigma=1, cov=0, n=100)
In [ ]:
plt.figure(figsize=(10,7))
plt.grid()
plt.xlabel(r'$x_1$') , plt.ylabel(r'$x_2$')
plt.scatter(datosNoCorrel[:,0], datosNoCorrel[:,1], s=14**2, alpha=0.75)
In [ ]:
datosCorrel = gaussiana2D(prom=50, sigma=1, cov=0.9, n=100)
In [ ]:
plt.figure(figsize=(10,7))
plt.grid()
plt.xlabel(r'$x_1$') , plt.ylabel(r'$x_2$')
plt.scatter(datosCorrel[:,0], datosCorrel[:,1], s=14**2, alpha=0.75)
In [ ]:
num_ptos = len(datosCorrel)
prom = datosCorrel.sum(0)/num_ptos # sum(0) suma sobre las columnas
print prom
In [ ]:
datosCorrelPromCero = datosCorrel - prom
print datosCorrel[0,:]
print datosCorrelPromCero[0,:]
Si $\scriptsize \mathbf{X} \in \mathbb{R}^{n \times d}$ es nuestra matriz de datos centrados en el origen, la matriz de covarianza está dada por:
$$ \mathbf{C}_{\mathbf X} = \frac{1}{n} \mathbf{X}^\top \mathbf{X} \,.$$Para que lo que estamos haciendo pueda ser paralelizable (en un momento explicaré por qué queremos ésto), en lugar de calcular directamente el producto de $\scriptsize \mathbf{X}^\top \mathbf{X}$ vamos a calcular el producto de estas dos matrices como la suma de productos exteriores.
In [ ]:
covCorrel = sum( map(lambda pto: np.outer(pto, pto), datosCorrelPromCero) ) / float(num_ptos)
print covCorrel
In [ ]:
def covarianza(datos):
## Llenar
Podemos usar la matriz de covarianza para encontrar las direcciones de máxima variación. Ésto se hace descomponiendo a la matriz para encontrar sus eigenvalores y eigenvectores. Una matriz en $\scriptsize \mathbb{R}^{d \times d}$, tiene $d$ eigenvalores y eigenvectore (no necesariamente distintos).
Recordamos que $\scriptsize \lambda$ es un eigenvalor de $\scriptsize \mathbf{X} \in \mathbb{R}^{d \times d}$ si existe $\scriptsize \mathbf{v}_{\lambda}$ tal que:
$$ \mathbf{X} \, \mathbf{v}_{\lambda} = \lambda \, \mathbf{v}_{\lambda}$$Para cada eigenvalor $\scriptsize \lambda$ existen infinidad de vectoress que satisfacen la ecuación de arriba, por lo que por conveción tomamos $\scriptsize \mathbf{v}_{\lambda}$ unitario.
Nos interesa quedarnos con las $k$ direcciones que concentran la mayor variación, esto es, queremos encontrar los $k$ eigenvectores asociados a los $k$ eigenvalores más grandes.
In [76]:
from numpy.linalg import eigh
eigVals, eigVecs = eigh(covCorrel)
print eigVals, '\n'
print eigVecs
In [ ]:
inds = np.argsort(eigVals)
top = eigVecs[:,inds[-1]]
print top
In [74]:
scoresCorrel = map(lambda pto: pto.dot(top), datosCorrel)
print scoresCorrel[0:3]
In [ ]:
def pca(datos, k=2):
cov = covarianza(datos)
eigVals, eigVecs = eigh(cov)
# Ordenar y sólo tomar k eigenVals y eigenVecs
# Calcular calificaciones
In [ ]:
def varianzaExplicada(datos, k=1):
components, scores, eigenvalues = pca(data,k)
cum = eigenvalues.cumsum()
return cum[k-1]/cum[-1]
In [ ]:
noCorrel1 = varianzaExplicada(datosNoCorrel, 1)
noCorrel2 = varianzaExplicada(datosNoCorrel, 2)
correl1 = varianzaExplicada(datosCorrel, 1)
correl2 = varianzaExplicada(datosCorrel, 2)
print 100*noCorrel1
print 100*noCorrel2, \n
print 100*correl1
print 100*correl2