Análisis de componentes principales con python

Esta entrada también se podría haber llamado:

'Reducción de la dimensión'

'Comprimiendo imágenes usando análisis de componentes principales y python'

Para la siguiente entrada se ha usado python 2.7.2, numpy 1.6.1, matplotlib 1.1.0 y sklearn 0.10

El análisis de componentes principales (PCA, por sus siglas en inglés, Principal Component Analysis) es una técnica que trata de reducir el número de dimensiones (número de variables) de un conjunto de datos intentando, a su vez, conservar la mayor cantidad de información. Es una técnica extremadamente útil como análisis exploratorio de datos (exploratory data analysis, pongo algunos términos en inglés porque a veces pueden resultar extraños en castellano), cuando se tiene demasiada información (muchas dimensiones, variables) y no se puede analizar correctamente la información. Se ha usado de forma exitosa para encontrar patrones, determinar 'outliers', compresión de imágenes,... Más adelante se explicará brevemente el proceso matemático que hay detrás.

Para el presente artículo vamos a usar el PCA para obtener imágenes comprimidas. Vamos a trabajar con una imagen monocromática extraída de la galería de AdeRussell en flickr. (con licencia cc 2.0) La imagen en cuestión es la siguiente:

Martin Luther King Memorial and the Washington Monument por AdeRussell
Martin Luther King Memorial and the Washington Monument por AdeRussell

Primero la guardaremos en nuestro equipo usando la biblioteca urllib2 disponible en la biblioteca estándar de vuestra instalación de python:

import urllib2
## Leemos la imagen desde la url
url = 'http://farm9.staticflickr.com/8312/8059125132_bed732fcc8_z.jpg'
kk = urllib2.urlopen(url).read()
## Guardamos la imagen en el directorio donde nos encontremos
## con el nombre 'king.jpg'
imagen = open('king.jpg', 'wb')
imagen.write(kk)
imagen.close()

La imagen está en formato jpg, la vamos a leer usando matplotlib y vamos a guardar la información en una matriz 2D (tenemos tres canales (r,g,b) que son iguales (imagen en escala de grises) por lo que solo vamos a usar uno de ellos).

import numpy as np
import matplotlib.pyplot as plt
## Leemos la imagen como un numpy array
kk = plt.imread('king.jpg')
## Si hacemos kk.shape vemos que existen
## tres canales en la imagen (r, g, b)
## Pero como es una imagen en escala de grises
## Los tres canales tienen la misma información
## por lo que nos podemos quedar con un solo canal
plt.subplot(221)
plt.title('canal 1')
plt.imshow(kk[:,:,0])
plt.subplot(222)
plt.title('canal 2')
plt.imshow(kk[:,:,1])
plt.subplot(223)
plt.title('canal 3')
plt.imshow(kk[:,:,2])
## Vemos que la imagen está rotada, hacemos uso de np.flipud
## http://docs.scipy.org/doc/numpy/reference/generated/numpy.flipud.html
plt.subplot(224)
plt.title('canal 1 rotado en BN')
plt.imshow(np.flipud(kk[:,:,0]), cmap=plt.cm.Greys_r)
plt.show()
## Finalmente, nos quedamos con una única dimensión
## Los tres canales rgb son iguales (escala de grises)
matriz = np.flipud(kk[:,:,0])

Bueno, ahora ya tenemos nuestra imagen como la queremos para poder empezar a trabajar con ella.

¿Cómo se obtienen las componentes principales?

Continue reading

La importancia de inspeccionar los datos

Muchas veces, cuando analizo datos, se me cuelan datos erróneos que me llevan a dolores de cabeza posteriores. Para evitar o minimizar esto, una buena práctica consiste en realizar una representación previa de esos datos para ver si los resultados posteriores son coherentes. Un buen ejercicio para entender la importancia de esto es echarle un ojo al cuartero de Anscombe. En el siguiente script se representa el cuarteto de Anscombe y sus principales estadísticos:

import numpy as np
import matplotlib.pyplot as plt
plt.ion()
ans_x_I = np.array([10.0, 8.0, 13.0, 9.0, 11.0,
                    14.0, 6.0, 4.0, 12.0, 7.0, 5.0])
ans_y_I = np.array([8.04, 6.95, 7.58, 8.81, 8.33,
                    9.96, 7.24, 4.26, 10.84, 4.82, 5.68])
ans_x_II = np.array([10.0, 8.0, 13.0, 9.0, 11.0,
                     14.0, 6.0, 4.0, 12.0, 7.0, 5.0])
ans_y_II = np.array([9.14, 8.14, 8.74, 8.77, 9.26,
                     8.10, 6.13, 3.10, 9.13, 7.26, 4.74])
ans_x_III = np.array([10.0, 8.0, 13.0, 9.0, 11.0,
                      14.0, 6.0, 4.0, 12.0, 7.0, 5.0])
ans_y_III = np.array([7.46, 6.77, 12.74, 7.11, 7.81,
                      8.84, 6.08, 5.39, 8.15, 6.42, 5.73])
ans_x_IV = np.array([8.0, 8.0, 8.0, 8.0, 8.0,
                     8.0, 8.0, 19.0, 8.0, 8.0, 8.0])
ans_y_IV = np.array([6.58, 5.76, 7.71, 8.84, 8.47,
                     7.04, 5.25, 12.50, 5.56, 7.91, 6.89])
x = [ans_x_I, ans_x_II, ans_x_III, ans_x_IV]
y = [ans_y_I, ans_y_II, ans_y_III, ans_y_IV]
def pinta(x, y, grupo, pos):
    plt.subplot(pos)
    minimo = np.min([x, y])
    maximo = np.max([x, y])
    plt.xlim(minimo - 1, maximo + 1)
    plt.ylim(minimo - 1, maximo + 1)
    plt.plot(x[grupo], y[grupo], 'yo')
    a1, a0 = np.polyfit(x[grupo], y[grupo], 1)
    plt.plot(x[grupo], a0 + a1 * x[grupo], 'b')
    plt.text(minimo, maximo - 0.5, 'y=%5.3f+%5.3fx' % (a0, a1))
    plt.text(minimo, maximo - 2.0, 'r=%5.3f' % np.corrcoef(x[grupo], y[grupo])[0,1])
    plt.text(minimo, maximo - 3.5, 'x media=%5.3f' % np.mean(x[grupo]))
    plt.text(minimo, maximo - 5, 'y media=%5.3f' % np.mean(y[grupo]))
    plt.text(minimo, maximo - 6.5, 'x std=%5.3f' % np.var(x[grupo]))
    plt.text(minimo, maximo - 8.0, 'y std=%5.3f' % np.var(y[grupo]))
pinta(x, y, 0, 221)
pinta(x, y, 1, 222)
pinta(x, y, 2, 223)
pinta(x, y, 3, 224)

Cuyo resultado sería el siguiente:

Se observa que los estadísticos representados son similares. Mismo ajuste lineal, coeficiente de correlación, media de los valores de x, media de los valores de y, varianza de los valores de x y varianza de los valores de y. Sin embargo, viendo los datos representados vemos que en el segundo no deberíamos usar un ajuste lineal para relacionar los valores de x y de y, en el tercer y cuarto grupos vemos que tenemos un dato que podría representar un dato erróneo (una mala medida, por ejemplo) y que en el tercer caso, al no haberlo filtrado correctamente, obtenemos un ajuste que nos llevaría a errores apreciables en el ajuste de los datos y en el cuarto caso nos llevaría a errores muy grandes en el ajuste de los datos mostrando una relación lineal que no parece existir.

Gracias a este plagio que acabo de hacer del artículo de la Wikipedia (traducido a python :-)) vemos la importancia de tratar nuestros datos con mimo, de realizar exploraciones iterativas y de refinar el conjunto final de datos 'limpios de impurezas' para hacer un análisis correcto de la 'realidad'. También vemos la importancia de compartir los datos base de los experimentos así como el código usado para el análisis de los mismos para que otros puedan reproducir los resultados y comprobar la veracidad de los mismos o poder corregir el análisis en caso necesario.

Finalmente os dejo un delirante vídeo de una charla TED de Hans Rosling donde se muestra la necesidad de mostrar de forma creativa los datos para poder llegar a un mejor entendimiento del universo que nos rodea.

[ted id=92 lang=es width=560 height=315]

P.D.: Sirva esta entrada para volver a homenajear a John Hunter y darle las gracias por permitirnos, mediante matplotlib, separarnos de los árboles para poder ver el bosque.

P.D.2: ¡¡¡¡¡No dejéis de ver el vídeo!!!!!

Estadística en Python con SciPy (I)

Introducción

Hoy vamos a ver cómo trabajar con variable aleatoria con el módulo stats de la biblioteca Scipy. Scipy viene con numerosas distribuciones de probabilidad, tanto discretas como continuas, y además pone a nuestra disposición herramientas para crear nuestras propias distribuciones y multitud de herramientas para hacer cálculos estadísticos. En esta primera parte nos centraremos en cómo manejar esas distribuciones y sus funciones de distribución, cómo representarlas con matplotlib y cómo definir nuevas distribuciones.

En esta entrada se ha usado python 2.7.3, numpy 1.6.1, matplotlib 1.1.0 y scipy 0.10.1.

Continue reading