import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

"""
Explanation: Principal Component Analysis
De Wikipedia:
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.
End of explanation
"""


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)

datosNoCorrel = gaussiana2D(prom=50, sigma=1, cov=0, n=100)

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)

datosCorrel = gaussiana2D(prom=50, sigma=1, cov=0.9, n=100)

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)

"""
Explanation: Vamos a generar datos:
End of explanation
"""


num_ptos = len(datosCorrel)
prom = datosCorrel.sum(0)/num_ptos # sum(0) suma sobre las columnas

print prom

datosCorrelPromCero = datosCorrel - prom

print datosCorrel[0,:]
print datosCorrelPromCero[0,:]

"""
Explanation: Centramos los datos:
End of explanation
"""


covCorrel = sum( map(lambda pto: np.outer(pto, pto), datosCorrelPromCero) ) / float(num_ptos)

print covCorrel

"""
Explanation: Calculamos la matriz de covarianza:
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.
End of explanation
"""


def covarianza(datos):
    ## Llenar

"""
Explanation: Ejercicio:
Escribir una función para calcular la covarianza de un conjunto de datos. Recibe como argumento un arreglo de puntos y regresa la matriz de covarianza.
End of explanation
"""


from numpy.linalg import eigh

eigVals, eigVecs = eigh(covCorrel)

print eigVals, '\n'
print eigVecs

"""
Explanation: Descomposición en eigenvalores:
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.
End of explanation
"""


inds = np.argsort(eigVals)
top = eigVecs[:,inds[-1]]

print top

"""
Explanation: Ordenamos a los eigenvalores:
End of explanation
"""


scoresCorrel = map(lambda pto: pto.dot(top), datosCorrel)

print scoresCorrel[0:3]

"""
Explanation: Calculamos  las "calificaciones" (scores) de los datos originales:
End of explanation
"""


def pca(datos, k=2):
    cov = covarianza(datos)
    
    eigVals, eigVecs = eigh(cov)
    # Ordenar y sólo tomar k eigenVals y eigenVecs
    
    # Calcular calificaciones

"""
Explanation: Ejercicio:
Escribir la función que hace el análisis de PCA entero. Recibe como argumento los datos y k el número de componentes principales que uno va a conservar. Regresa las k componentes principales (los eigenvectores), las calificaciones de los datos originales, y los k eigenvalores principales.
End of explanation
"""


def varianzaExplicada(datos, k=1):
    components, scores, eigenvalues = pca(data,k)
    
    cum = eigenvalues.cumsum()
    
    return cum[k-1]/cum[-1]

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

"""
Explanation: ¿Qué porcentaje de la varianza se explica con k componentes?
End of explanation
"""


