Importamos las librerias necesarias:
%matplotlib inline
# importamos las librerías básicas
import random
import numpy as np
import matplotlib.pyplot as plt
# importamos librerías y class para dibujar los vectores que definen las dimensiones
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
self._verts3d = xs, ys, zs
def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
FancyArrowPatch.draw(self, renderer)
Creamos un conjunto de muestras con las que trabajar en el siguiente ejemplo. Distribución normal de 100 muestras con 3 características cada una:
# centro de la distribución
mean = np.array([1,1,1.5])
# matriz de covarianza de la distribución
cov = np.array([[1.2,0,0],[0,.9,0],[0,0,1.8]])
# número de muestras
size = 100
# creamos la distribución
X = np.random.multivariate_normal(mean, cov, size)
Dibujamos el conjunto de muestras creado:
fig = plt.figure()
ax = fig.gca(projection="3d")
ax.plot(X[:,0], X[:,1], X[:,2], "ro")
ax.set_title("Muestras X")
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_zlabel("x3")
plt.show()
Tratamos los datos, normalizamos y escalamos:
# normalizamos los datos
X -= np.mean(X, 0)
# escalamos las características
X /= np.std(X, 0)
# escalamos todos los datos a 1
X /= np.abs(X).max()
Dibujamos el nuevo conjunto de datos optimizado:
fig = plt.figure()
ax = fig.gca(projection="3d")
ax.plot(X[:,0], X[:,1], X[:,2], "co")
ax.set_title("Muestras X escaladas")
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_zlabel("x3")
plt.show()
Aplicamos el algoritmo PCA (Principal Component Analysis).
Primero calculamos e imprimimos la matriz de covarianza de los datos optimizados:
# centro de la distribución
mean_f1 = np.mean(X[:,0])
mean_f2 = np.mean(X[:,1])
mean_f3 = np.mean(X[:,2])
# matriz covarianza
C = np.cov(X.T)
print "Matriz de covarianza de X:"
print C
Matriz de covarianza de X: [[ 0.14941101 0.01186524 0.00770462] [ 0.01186524 0.14941101 -0.01495807] [ 0.00770462 -0.01495807 0.14941101]]
Calculamos autovectores y autovalores de la matriz de covarianza con la función Descomposición en valores singulares svd() (autovectores ordenados según autovalores de mayor a menor):
U, s, V = np.linalg.svd(C)
print "Autovalores:"
print s
print "Autovectores:"
print U
Autovalores: [ 0.16526682 0.15682372 0.12614251] Autovectores: [[-0.28938195 0.80722322 -0.51444024] [-0.76247841 0.13052089 0.63371206] [ 0.57869229 0.5756344 0.57771989]]
Dibujamos los autovectores sobre los datos:
fig = plt.figure()
ax = fig.gca(projection="3d")
ax.plot(X[:,0], X[:,1], X[:,2], "co")
ax.plot([0], [0], [0], "ko")
for v in U.T:
a = Arrow3D([0, v[0]], [0, v[1]], [0, v[2]], mutation_scale=20, arrowstyle="-|>", color="k")
ax.add_artist(a)
ax.set_title("Muestras X y sus autovectores")
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_zlabel("x3")
plt.show()
Calculamos las nuevas 'k' dimensiones según la varianza que queramos mantener de los datos originales. Al ser un ejemplo, con el objetivo de mostrar la reducción de dimensionalidad, se mantiene una varianza muy baja:
def dim_reduction(s, var):
for k in range(1,1000):
variance_retaines = sum(s[0:k]) / sum(s)
if variance_retaines >= var:
return k
k = dim_reduction(s, .69)
print "Reducción de dimensiones. De 3D a " + str(k) + "D"
Reducción de dimensiones. De 3D a 2D
Dibujamos la proyección de los datos en el nuevo espacio 2D:
U_reduced = U[:,:k]
Z = U_reduced.T.dot(X.T)
Z = -Z.T
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(Z[:,0], Z[:,1], "go")
ax.set_title("Datos en el nuevo espacio")
ax.set_xlabel("z1")
ax.set_ylabel("z2")
plt.show()