Regresión en Python

La regresión.

En varias entradas de Series de Tiempo y de Machine Learning en R comenté sobre la regresión lineal. El método es bien conocido en varios campos y quizás es el más natural de aprender, ya que la idea de trazar la recta que mejor describe la relación entre puntos de un plano parece familiar e intuitivo.

Ejemplo_opt_regresion

El problema con el método o el modo en que se enseña pude variar, como el método en ocasiones se subvalora debido a que ese considera en general que es solo “agregar línea de tendencia” a los datos, pero se deja de lado la revisión de los estadísticos y de revisar el comportamiento de los errores de la estimación.

En la entrada de manera breve comento un pequeño recorrido sobre las técnicas de regresión clásica y las penalizadas. Las principales referencias sobre las cuales comento los métodos son [1,2] y algunos artículos o reportes que dejo en las referencias. Lo que ya no comento es la parte estadística; es decir, los detalles sobre las pruebas de hipótesis y las potencias de las pruebas las comenté en otras entradas. En esta trato de solo mencionar lo mínimo sobre esos detalles, pero enfocarme más en otros usos y las variaciones a la técnica tradicional.

El modelo en general.

La idea de la regresión es tener un conjunto de datos de entrada y explicar o describir las salidas como una combinación lineal de los datos de entrada.

Es decir, se tienen los datos: Datos_entrada

y el modelo que se pretende describir es:

Modelo_Lineal

En la ecuación anterior se interpreta como la salida del modelo f(X) sea lo más parecida a Y, que son los valores reales a estimar.

Para el caso de los mínimos cuadrados se tiene que la función a minimizar es:

LQ

LQ2

Lo que se hace es sustituir la expresión que de f(x) que depende de las β’s para minimizar la ecuación, lo que se hace es tratar de encontrar los mejores β’s para describir Y.

En los métodos modernos lo que se penalizar la estimación de los β’s, debido a que al estimar la regresión por mínimos cuadrados sobre un subconjunto de entrenamiento ( train set) de los datos, es decir algunos de los x’s, permite obtener los parámetros  β’s sobre estimados  y al momento de usar el modelo para el conjunto de datos de prueba (test set) se obtienen tasas de errores mayores.

Los métodos con penalización usuales son Ridge y Lasso, los dos agregan un parámetro λ que implica el problema de tener un método para definir el mejor valor para el modelo. Explico con ejemplos que se hace para determinar el parámetro.

Para el modelo Ridge a la función de errores es como la de los mínimos cuadrados pero agregando un término más. La ecuación es la siguiente:

RR

Para el modelo Lasso es similar a la de Ridge, pero en la penalización en lugar de considerar los β’s  al cuadrado se considera su valor absoluto. La ecuación es la siguiente:

LR

Lamento poner las ecuaciones, pero creo que pese a que no se entienda la teoría de manera total , al observando las diferencias en las ecuaciones se puede apreciar la diferencia en los métodos. En la primera gráfica de la entrada no lo mencioné, pero si se presiona sobre la imagen y se observa con cuidado que tiene trazadas 3 rectas sobre los datos de color negro, estas tres rectas corresponden a los 3 métodos distintos.

Una visión o idea gráfica de lo que pasa con el modelo o técnica, es visto en la siguiente imagen. Al final el objetivo de la regresión lineal es encontrar el plano que mejor satisface ciertas condiciones para asegurar que predice o explica lo mejor posible los datos.

Regression_Lineal

Imagen tomada del libro de Trevor Hastie y Robert Tibshirani.

Ejemplo 1

Para visualizar la diferencia en las estimaciones haciendo uso de Python, sobre todo de la librería scikit-learn hago uso del código de Jaques Grobler, modifico parte de su código para agregar la estimación de la regresión Rigde y Lasso, además elijo una muestra mayor de datos para prueba para que la gráfica sea notoriamente más atractiva.

Lo que hace el código:

  1. Se toma una base de datos sobre diabetes la cual se puede cargar directamente en la librería scikit-learn.
  2. Se divide la muestra en datos de entrenamiento( train) y de prueba (test).
  3. Se estiman los 3 modelos de regresión y se estiman los estadísticos base de las estimaciones.
  4. Se hace la gráfica de los datos de test.
 # -*- coding: 850 -*-
print(__doc__)


# Code source: Jaques Grobler
# License: BSD 3 clause

#Se cargan los modulos

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

# Se cargan los datos
diabetes = datasets.load_diabetes()


# Se usa cierta muestra de los datos
diabetes_X = diabetes.data[:, np.newaxis]
diabetes_X_temp = diabetes_X[:, :, 2]

# Se separa la muestra de datos
diabetes_X_train = diabetes_X_temp[:-30]
diabetes_X_test = diabetes_X_temp[-30:]

# Se separa las variables explicadas
diabetes_y_train = diabetes.target[:-30]
diabetes_y_test = diabetes.target[-30:]

# Se crean los modelos lineales
regr = linear_model.LinearRegression()
regr2=linear_model.Lasso(alpha=.5)
regr3=linear_model.Ridge(alpha=.5)

#Se entrenan los modelos
regr.fit(diabetes_X_train, diabetes_y_train)
regr2.fit(diabetes_X_train, diabetes_y_train)
regr3.fit(diabetes_X_train, diabetes_y_train)


print u'Regresión Mínimos Cuadrados Ordinarios'
# Coeficiente
print'Coeficientes:',regr.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Lasso' 
# Coeficiente
print'Coeficientes:', regr2.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr2.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr2.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Ridge'
#Coeficiente
print'Coeficientes:', regr3.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr3.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr3.score(diabetes_X_test, diabetes_y_test))


# Plot 
plt.scatter(diabetes_X_test, diabetes_y_test, color='black')
plt.plot(diabetes_X_test, regr.predict(diabetes_X_test), color='blue',
 linewidth=3, label=u'Regresión MCO')
plt.plot(diabetes_X_test, regr2.predict(diabetes_X_test), color='yellow',
 linewidth=3, label=u'Regresión Lasso')
plt.plot(diabetes_X_test, regr3.predict(diabetes_X_test), color='green',
 linewidth=3, label=u'Regresión Ridge')
plt.title(u'Regresión Lineal por 3 metodos diferentes')
plt.legend()
plt.xticks(())
plt.yticks(())

plt.show()

Los resultados serían los siguientes:

Regresión Mínimos Cuadrados Ordinarios
Coeficientes: [ 941.43097333]
Residual sum of squares: 3035.06
Varianza explicada: 0.41

Regresión Lasso
Coeficientes: [ 720.32135439]
Residual sum of squares: 3265.36
Varianza explicada: 0.37

Regresión Ridge
Coeficientes: [ 612.64202818]
Residual sum of squares: 3458.17
Varianza explicada: 0.33

La gráfica se vería de la siguiente forma:

Tres_Métodos

El ejemplo anterior viendo la gráfica y los estadísticos tenemos como conclusión que le método MCO; mínimo cuadrados ordinarios, es el que mejor se ajusta a nuestros datos. En este caso uno se pregunta sobre la ventaja de usar otro método para la regresión lineal, si los resultados son favorables para la técnica más conocida.

Lo que no se menciona en el código, o lo que se queda sin explicar es la elección de λ para los modelos Ridge y Lasso. Lo que se hace para aplicar un buen modelo de estos dos métodos, es estrenar el modelo de modo tal que el valor de λ sea el “mejor”. Para esto se usa Cross-Validation, así que haciendo ligeros cambios al código anterior se puede mejorar la aplicación de los modelos lineales.

 # -*- coding: 850 -*-
print(__doc__)


# Code source: Jaques Grobler
# License: BSD 3 clause

#Se cargan los módulos
import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

# Se cargan los datos
diabetes = datasets.load_diabetes()


# Se eligen las muestras de datos
diabetes_X = diabetes.data[:, np.newaxis]
diabetes_X_temp = diabetes_X[:, :, 2]

# Se dividen los datos en training/testing 
diabetes_X_train = diabetes_X_temp[:-30]
diabetes_X_test = diabetes_X_temp[-30:]

# Se divide la variable explicada en  training/testing 
diabetes_y_train = diabetes.target[:-30]
diabetes_y_test = diabetes.target[-30:]

# Se crean los modelos lineales
regr = linear_model.LinearRegression()
#regr2=linear_model.LassoCV(alpha=)
regr3=linear_model.RidgeCV(alphas=[0.1,0.2,0.5,1.0,3.0,5.0,10.0])

# Se entrenan los modelos
regr.fit(diabetes_X_train, diabetes_y_train)
regr2=linear_model.LassoCV(cv=20).fit(diabetes_X_train, diabetes_y_train)
regr3.fit(diabetes_X_train, diabetes_y_train)


print u'Regresión Mínimos Cuadrados Ordinarios'
# Coeficiente
print'Coeficientes:',regr.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Lasso' 
# Coeficientes
print'Coeficientes:', regr2.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr2.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr2.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Ridge'
# Coeficiente
print'Coeficientes:', regr3.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr3.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr3.score(diabetes_X_test, diabetes_y_test))


# Plot 
plt.scatter(diabetes_X_test, diabetes_y_test, color='black')
plt.plot(diabetes_X_test, regr.predict(diabetes_X_test), color='blue',
 linewidth=3, label=u'Regresión MCO')
plt.plot(diabetes_X_test, regr2.predict(diabetes_X_test), color='yellow',
 linewidth=3, label=u'Regresión Lasso')
plt.plot(diabetes_X_test, regr3.predict(diabetes_X_test), color='green',
 linewidth=3, label=u'Regresión Ridge')
plt.title(u'Regresión Lineal por 3 metodos diferentes')
plt.legend()
plt.xticks(())
plt.yticks(())

plt.show()

Lo que se tienen como resultado es lo siguiente:

Regresión Mínimos Cuadrados Ordinarios
Coeficientes: [ 941.43097333]
Residual sum of squares: 3035.06
Varianza explicada: 0.41

Regresión Lasso
Coeficientes: [ 940.48954236]
Residual sum of squares: 3035.57
Varianza explicada: 0.41

Regresión Ridge
Coeficientes: [ 850.17738251]
Residual sum of squares: 3103.11
Varianza explicada: 0.40

La gráfica correspondiente es:

Tres_Métodos2

 

Lo que se aprecia es que los modelos mejoran considerablemente, si bien el modelo MCO sigue siendo uno de los mejores para este ejemplo, algo a considerar es que es una muestra pequeña y de pocas dimensiones. Así que lo interesante es aplicar la regresión a muestras mayores y de datos con una mayor cantidad de variables.

Ejemplo 2-Regresión Multiple

Para el siguiente ejemplo hago uso de una matriz de datos de tamaño 506 líneas de datos con 12 variables de entrada y una variable que se trata de estimar por métodos lineales. Los datos corresponden a los valores de las viviendas en los suburbios de Boston.

Tomo estos datos para mostrar como aplicar la regresión a datos con varia variables, para este caso no hago una gráfica pero se pueden revisar los datos o valores de los estadísticos para tener un modo de elegir entre modelos.

El código es muy similar al anterior y lo que se puede hacer como variación es aplicar alguna técnica de reducción de dimensiones para elegir las variables que “impactan realmente” en el modelo o jugar con ejecutar la regresión no con todas las variables, sino empezar con unas y poco a poco agregar nuevas para revisar como se comporta el modelo al estimar con solo unas cuantas variables. Algo que tampoco hago en este caso es explorar las variables antes de aplicar la regresión, pero la intención es ver que prácticamente es igual la ejecución de una regresión multi lineal.

#import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

#13 variables explicativas

# Se cargan los datos
Boston = datasets.load_boston()



# Se divide en training/testing sets
Boston_train = Boston.data[:405,]
Boston_test = Boston.data[405:,]


# La variable explicativa se divide en training/testing sets
Boston_y_train = Boston.target[:405]
Boston_y_test = Boston.target[405:]

regr = linear_model.LinearRegression()
regr2=linear_model.LassoCV(cv=20).fit(Boston_train, Boston_y_train)
regr3=linear_model.RidgeCV(alphas=[0.1,0.2,0.5,0.7,1.0,1.5,3.0,5.0,7.0,10.0])
#np.size(Boston_y_train = Boston.target[:405])
regr.fit(Boston_train, Boston_y_train)
regr3.fit(Boston_train,Boston_y_train)

print u'Regresión Mínimos Cuadrados Ordinarios'
#Coeficiente
print'Coeficientes:',regr.coef_
# MSE 
print("Residual sum of squares: %.2f"
 % np.mean((regr.predict(Boston_test) - Boston_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr.score(Boston_train, Boston_y_train))


print u'Regresión Ridge'
# Coeficientes
print'Coeficientes:', regr3.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr3.predict(Boston_test) - Boston_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr3.score(Boston_train, Boston_y_train))


print u'Regresión Lasso' 
# Coeficiente
print'Coeficientes:', regr2.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr2.predict(Boston_test) - Boston_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr2.score(Boston_train, Boston_y_train))

Los resultados del código son los siguientes resultados:

Regresión Mínimos Cuadrados Ordinarios
Coeficientes: [ -1.94651664e-01 4.40677436e-02 5.21447706e-02 1.88823450e+00
 -1.49475195e+01 4.76119492e+00 2.62339333e-03 -1.30091291e+00
 4.60230476e-01 -1.55731325e-02 -8.11248033e-01 -2.18154708e-03
 -5.31513940e-01]
Residual sum of squares: 33.58
Varianza explicada: 0.74

Regresión Ridge
Coeficientes: [ -1.93961816e-01 4.42902027e-02 4.66864588e-02 1.87777854e+00
 -1.37158180e+01 4.76760163e+00 1.65803755e-03 -1.28477100e+00
 4.57255634e-01 -1.57415769e-02 -7.97849246e-01 -1.78511420e-03
 -5.33113498e-01]
Residual sum of squares: 33.15
Varianza explicada: 0.74

Regresión Lasso
Coeficientes: [-0.13537099 0.04758408 -0. 0. -0. 3.26851037
 0.01014632 -0.9131925 0.39587576 -0.01796467 -0.65864709 0.00511243
 -0.67828222]
Residual sum of squares: 23.48
Varianza explicada: 0.71

Si se tienen cuidado se puede observar que los dos primeros modelos muestran mucha “cercanía”, los dos tienen estimaciones casi similares y estadísticos muy parecidos. El modelo Lasso al contrario, parece considerablemente diferente y con valores considerablemente diferentes. El punto clave es cómo elegir entre los modelos. En esta parte sería adecuado hacer uso de algunos indicadores sobre la calidad del modelo como AIC o BIC, etc.

Ejemplo 3. Usando la regresión de otro modo.

El siguiente ejemplo es un ejemplo que considero bastante bonito y muestra el uso de Lasso  y Ridge que en general no se piensa, respecto el código y el reconocimiento a Emmanuelle Gouillart, sobre la cual recomiendo revisar su tutorial sobre procesamiento de imágenes. No hago cambio alguno al código, lo que hago es explicar que hace en cada etapa para clarificar y valorar las regresiones con penalización como herramientas de reconstrucción de imágenes.

El código:

print(__doc__)

# Author: Emmanuelle Gouillart <emmanuelle.gouillart@nsup.org>
# License: BSD 3 clause

import numpy as np
from scipy import sparse
from scipy import ndimage
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt


def _weights(x, dx=1, orig=0):
 x = np.ravel(x)
 floor_x = np.floor((x - orig) / dx)
 alpha = (x - orig - floor_x * dx) / dx
 return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))


def _generate_center_coordinates(l_x):
 l_x = float(l_x)
 X, Y = np.mgrid[:l_x, :l_x]
 center = l_x / 2.
 X += 0.5 - center
 Y += 0.5 - center
 return X, Y


def build_projection_operator(l_x, n_dir):
 """ Compute the tomography design matrix.

 Parameters
 ----------

 l_x : int
 linear size of image array

 n_dir : int
 number of angles at which projections are acquired.

 Returns
 -------
 p : sparse matrix of shape (n_dir l_x, l_x**2)
 """
 X, Y = _generate_center_coordinates(l_x)
 angles = np.linspace(0, np.pi, n_dir, endpoint=False)
 data_inds, weights, camera_inds = [], [], []
 data_unravel_indices = np.arange(l_x ** 2)
 data_unravel_indices = np.hstack((data_unravel_indices,
 data_unravel_indices))
 for i, angle in enumerate(angles):
 Xrot = np.cos(angle) * X - np.sin(angle) * Y
 inds, w = _weights(Xrot, dx=1, orig=X.min())
 mask = np.logical_and(inds >= 0, inds < l_x)
 weights += list(w[mask])
 camera_inds += list(inds[mask] + i * l_x)
 data_inds += list(data_unravel_indices[mask])
 proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
 return proj_operator


def generate_synthetic_data():
 """ Synthetic binary data """
 rs = np.random.RandomState(0)
 n_pts = 36.
 x, y = np.ogrid[0:l, 0:l]
 mask_outer = (x - l / 2) ** 2 + (y - l / 2) ** 2 < (l / 2) ** 2
 mask = np.zeros((l, l))
 points = l * rs.rand(2, n_pts)
 mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
 mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)
 res = np.logical_and(mask > mask.mean(), mask_outer)
 return res - ndimage.binary_erosion(res)


# Generate synthetic images, and projections
l = 128
proj_operator = build_projection_operator(l, l / 7.)
data = generate_synthetic_data()
proj = proj_operator * data.ravel()[:, np.newaxis]
proj += 0.15 * np.random.randn(*proj.shape)

# Reconstruction with L2 (Ridge) penalization
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)

# Reconstruction with L1 (Lasso) penalization
# the best value of alpha was determined using cross validation
# with LassoCV
rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)

plt.figure(figsize=(8, 3.3))
plt.subplot(131)
plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
plt.axis('off')
plt.title('original image')
plt.subplot(132)
plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L2 penalization')
plt.axis('off')
plt.subplot(133)
plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L1 penalization')
plt.axis('off')

plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
 right=1)

plt.show()

La imagen que se obtiene es la siguiente:

Img1_Regression

Primero explicando el resultado lo que se tienen es una imagen que se crea y se toma como imagen original, la imagen del centro y la derecha son aproximaciones lineales a la original. Lo interesante es que las dos aproximaciones son buenas, pero la generada por Lasso es mucho mejor que la generada por Ridge.

¿Cómo hace esto el código? En general la idea de aplicar una regresión es para predecir datos y se piensa en la imagen de la recta o el plano que se aproxima a los datos que se tratan de modelar.

En esta caso también se hace una aproximación lineal, lo cual es generado por las funciones:

1.-generate_synthetic_data()
2.-build_projection_operator()

Entonces en resumen, lo primero que se hace es cargan las librerías requeridas, lo cual es normal en todos los códigos. Luego se definen funciones que hacen uso de las funciones base de numpy para procesar y construir matrices. Explotando esto, lo que se construyen las dos funciones mencionadas anteriormente.

Así la muestra de datos queda considerada dentro de la variable proj sobre la cual se aplica el método np.ravel(), con lo cual se aplican las dos regresiones penalizadas para  hacer la aproximación lineal con la imagen original que se gurda en la variable data.

Lamento no explicar todo a mucho detalle, pero creo que los comentario de Emmanuelle son claros y los pocos conceptos que están ocultos en el código son claros con un poco de conocimiento en álgebra lineal.

Espero que aprecien el uso de los dos métodos de regresión lineal penalizados siendo usados en una aproximación lineal a un conjunto de datos diferente a lo que suele uno ver en los libros de texto.

Ejemplo 4.-Otro ejemplo con ciertos datos.

En el libro de Willi Richert y Luis Pedro Coelho, dedican dos capítulos a mostrar el uso de la regresión lineal como una herramienta para  predecir el  rating  en un sistema de recomendación. Recomiendo leer los capítulos 7 y 8 para que sea claro como se hace uso de la regresión para éste tipo de sistemas, referencia [3]. Como ejemplo de como usar la regresión es bastante bueno,  pero algo largo y no creo bueno compartirlo en el blog ya que son varios script para mostrar como funcionaria el sistema. Por ello recomiendo comprar el texto y revisar el código que comparten en su repositorio en Github para replicar y estudiar dicho ejemplo.

Un punto clave para los sistemas de recomendación es considerar el caso cuando se tienen una cantidad pequeña de usuario comparada con una gran cantidad de productos. Ejemplo, suponiendo que se va el supermercado, se tienen una cantidad de quizás 5 mil artículo o productos y se reciben en dicha tienda solo mil clientes. Construyendo la matriz con mil filas y 5 mil columnas, se tienen un problema del tipo Ρ»Ν.

Este tipo de problemas suelen aparecer mucho cuando uno analiza datos de los empleados de una empresa, sobre productos o cierto tipo de inventarios o matrices de datos que tienen que ver con los procesos de algún empresa.

En caso de considerar la regresión como un buen modelo aplicar a este tipo de datos, el problema es que se tiene una matriz de datos que no permitiría hacer una división de entre train/ test para el algoritmo, entonces como alternativa lo que se hace es emplear una técnica de Cross-Validation la cual permite elegir los parámetros de modelo con “mejor calidad”.

La técnica es llamada k-fold cross-validation, lo que se hace es que se divide la muestra en k particiones y se elige una de esas como el conjunto de test y las k-1 como train. Se repite el proceso k veces y al final se eligen las k-medias de los coeficientes como los valores de los coeficientes del modelo. Dejo en las referencias los tutoriales de Andrew Moore, en el cual pueden ver en las presentaciones como funciona de manera general el método.

CV-kfold

Ejemplo gráfico de k-fold Croos-Validation

Sobre el código, entonces como ejemplo hago uso del código de Willi Richert y Luis Pedro Coelho con una ligera modificación para considerar el uso de Lasso y comparar los resultados.

#Módulos
import numpy as np
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.datasets import load_svmlight_file
from sklearn.linear_model import LinearRegression, LassoCV
from sklearn.cross_validation import KFold

#Cargar los datos
data, target = load_svmlight_file('data/E2006.train')

lr = LinearRegression()

# Scores:

#Estimación del modelo lineal
lr.fit(data, target)

#Estimación del modelo Lasso
lr2=LassoCV(cv=5).fit(data,target)

#Estimaciones o predicciones
pred = lr.predict(data)
pred2 = lr2.predict(data)


print('RMSE sobre los datos de entrenamiento, {:.2}'.format(np.sqrt(mean_squared_error(target, pred))))

print('R2 sobre los datos de entrenamiento, {:.2}'.format(r2_score(target, pred)))
print('')

pred = np.zeros_like(target)
kf = KFold(len(target), n_folds=5)
for train, test in kf:
 lr.fit(data[train], target[train])
 pred[test] = lr.predict(data[test])

print('RMSE sobre los datos de prueba (5-fold), {:.2}'.format(np.sqrt(mean_squared_error(target, pred))))
print('R2 sobre los datos de prueba (5-fold), {:.2}'.format(r2_score(target, pred)))

#Estimaciones para Lasso
print('RMSE en Lasso(5-fold), {:.2}'.format(np.sqrt(mean_squared_error(target, pred2))))
print('R2 en Lasso(5-fold), {:.2}'.format(r2_score(target, pred2)))

El código regresa cierta información, al cuales son los estadísticos básicos para saber como se comporta el modelo.

RMSE sobre los datos de entrenamiento, 0.0012
R2 sobre los datos de entrenamiento, 1.0

RMSE sobre los datos de prueba(5-fold), 0.59
R2 sobre los datos de prueba(5-fold), 0.12
RMSE en Lasso(5-fold), 0.37
R2 en Lasso(5-fold), 0.66

Si se revisan los primeros datos causa confusión, por que parecen que el primer modelo es “perfecto”, pero eso es una mala señal sobre el modelo. Para más detalles dejo algunas referencias sobre sobre los estadísticos y su valoración.

Entonces con este último ejemplo la intención es mostrar que cuando se tienen condiciones del tipo Ρ»Ν  lo adecuado es usar Cross-Validation para este tipo de muestras al momento de aplicar la regresión.

Varios errores que cometemos al aplicar una regresión lineal.

En un post de Cheng-Tau Chu listó 7 errores que pasan al aplicar algoritmos de Machine Learning. Eligiendo los errores que se presentan en la regresión lineal en general, la lista sería algo así:

  • No revisar la correlación de los errores
  • No tener varianza constante  en los errores.
  • No considerar la importancia de los outliers.
  • No cuidar los problemas de la multicolinealidad o colinealidad en los datos.
  • Considerar un modelo lineal cuando las interacciones de los datos son no-lineales.

Hablar de cada punto requiere dar un ejemplo de donde falla el modelo lineal y como abordar dicho problema. Recomiendo la referencia [2] para revisar algunos detalles.

Conclusión: Espero que los ejemplos den una idea general sobre como se usa la regresión lineal. Hay otras variaciones de la regresión que son LAR y Elastic Net, sobre las dos hay mucho material y sobre todo que son técnicas más nuevas y han estado en recientes cambios con nuevas investigaciones.

Referencias:

  1. https://mitpress.mit.edu/books/machine-learning-2
  2. http://statweb.stanford.edu/~tibs/ElemStatLearn/
  3. http://www.amazon.es/Building-Machine-Learning-Systems-Python/dp/1782161406
  4. http://www.springer.com/us/book/9780387310732

 

Anuncios

Suppor Vector Machine ( Máquina de Soporte de Vectores) en Python

Sobre SVM

Las técnicas de Máquina de Soporte de Vectores (SVM) por alguna extraña razón es uno de los algoritmos que más interés despiertan y sobre todo el que creo que genera más sobre estimación. Es posible que no he usado suficiente la familia de algoritmos como para reconocer tanto su valor.

Hace no más de un año un amigo en una propuesta sugería usar esta familia de algoritmos sin revisar si realmente convenía hacer uso de ellos para el tipo de datos que se analizarían. La justificación de este amigo era que SVM eran los “mejores” algoritmos, lo cual es falso, no se puede calificar a una familia de algoritmos como los mejores, se requiere probar varios contra los datos para determinar cuales funcionan mejor con ellos.

Dos ejemplos que he compartido donde involucro SVM se pueden revisar en las entradas:

1.-Comparacion-entre-svm-naive-bayes-arboles-de-decision-y-metodos-lineales.

2.-Maquina-de-soporte-vectorial-svm-sopport-vector-machine

En las dos entradas uso SVM para clasificar, si bien sus resultados son buenos en las muestras de datos no resulta ser el mejor método al compararlo con otra familia de algoritmos. Lo que recomiendo es probar con más de una familia de algoritmos para analizar alguno de los problemas y elegir más de un estadístico o indicador que puedan servir para comparar la eficiencia de los algoritmos.

Si se tiene interés en conocer más detalles sobre los algoritmos recomiendo consultar la referencia [1,2,3,4], en orden con respecto a la dificultad del contenido y de las explicaciones. Los primeros 3 ejemplos se pueden encontrar en la documentación de la librería scikit-learn.

La relación entre las técnicas o algoritmos es la siguiente; el primero es Margen Maximal de Clasificación (MMC), su generalización es Soporte de Vectores para Clasificación (SVC) y esta a su vez se generaliza a SVM.

El primer ejemplo con datos simulados, son la imagen clásica para explicar como funciona SVM, pero lo que se muestra en sí es el caso MMC.

El código es el siguiente:

#MMC 
print(__doc__)

#Librerías requeridas

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm

# Se crean los datos
np.random.seed(0)
X = np.r_[np.random.randn(70, 2) - [2, 2], np.random.randn(70, 2) + [2, 2]]
Y = [0] * 70 + [1] * 70

#Se estima el modelo
clf = svm.SVC(kernel='linear')
clf.fit(X, Y)

# Se construye la recta que separa las clases
w = clf.coef_[0]
a = -w[0] / w[1]
xx = np.linspace(-5, 5)
yy = a * xx - (clf.intercept_[0]) / w[1]


# support vectors
b = clf.support_vectors_[0]
yy_down = a * xx + (b[1] - a * b[0])
b = clf.support_vectors_[-1]
yy_up = a * xx + (b[1] - a * b[0])

#Gráfica
plt.plot(xx, yy, 'k-')
plt.plot(xx, yy_down, 'k--')
plt.plot(xx, yy_up, 'k--')

plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1],
 s=80, facecolors='none')
plt.scatter(X[:, 0], X[:, 1], c=Y, cmap=plt.cm.Paired)

plt.title('Recta separadora')
plt.axis('tight')
plt.show()

La gráfica es la siguiente:

MMC

A lo que se refiere el algoritmo con el nombre de Margen Maximal, es al margen que se obtiene entre las dos clases de datos, siendo el margen la distancia entre las líneas punteadas.

Esto en teoría es una optimización sobre el margen y la generalización se da cuando no se puede obtener una línea recta que no intercepte o que separe de manera clara las muestras de datos. El ejemplo clásico se hace con los datos Iris, los cuales permiten mostrar como se generaliza el uso de los margenes que separan los datos, cuando no son líneas rectas.

El código muestra el ejemplo de como se separan de manera no lineal las fronteras.

#Código

print(__doc__)

#Módulos
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets

#Se importan los datos
iris = datasets.load_iris()
X = iris.data[:, :2] 
y = iris.target

h = .02

#Parámetro de regularización
C = 1.0 

svc = svm.SVC(kernel='linear', C=C).fit(X, y)
rbf_svc = svm.SVC(kernel='rbf', gamma=0.7, C=C).fit(X, y)
poly_svc = svm.SVC(kernel='poly', degree=3, C=C).fit(X, y)
lin_svc = svm.LinearSVC(C=C).fit(X, y)

#Probé pero no me gustó el resultado con sigmoid
#sigmoid_svc=svm.SVC(kernel='sigmoid').fit(X,y)

#Se crean los marcos para las gráficas

x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h),np.arange(y_min, y_max, h))

#Títulos para las gráficas
titles = ['SVC con kernel lineal',
 'Lineal SVC',
 'SVC con RBF kernel',
 'SVC con polinomio(grado 3) kernel']


for i, clf in enumerate((svc,lin_svc, rbf_svc, poly_svc)):
 # Se grafican las fronteras 
 plt.subplot(2, 2, i + 1)
 plt.subplots_adjust(wspace=0.4, hspace=0.4)

 Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

 #Color en las gráficas
 Z = Z.reshape(xx.shape)
 plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.8)

 #Puntos de entrenamiento
 plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)
 plt.xlabel('Longitud Sepal')
 plt.ylabel('Peso Sepal')
 plt.xlim(xx.min(), xx.max())
 plt.ylim(yy.min(), yy.max())
 plt.xticks(())
 plt.yticks(())
 plt.title(titles[i])

plt.show()

La gráfica que se obtiene es la siguiente:

SVC_Iris

El otro uso es para estimar la regresión y un ejemplo visual que se aprecia rápido es en la regresión de datos en dos dimensiones, es decir; en el plano.

El código es el siguiente:

#Generación de datos

print(__doc__)
import numpy as np
from sklearn.svm import SVR
import matplotlib.pyplot as plt

#Generación de datos
X = np.sort(5 * np.random.rand(100, 1), axis=0)
y = np.sin(X).ravel()
y[::5] += 3 * (0.5 - np.random.rand(20))

#Modelos regresión fit

svr_rbf = SVR(kernel='rbf', C=1e3, gamma=0.1)
svr_lin = SVR(kernel='linear', C=1e3)
svr_poly = SVR(kernel='poly', C=1e3, degree=2)
y_rbf = svr_rbf.fit(X, y).predict(X)
y_lin = svr_lin.fit(X, y).predict(X)
y_poly = svr_poly.fit(X, y).predict(X)

#Mirando los resultados

plt.scatter(X, y, c='k', label='datos')
plt.hold('on')
plt.plot(X, y_rbf, c='g', label='Modelo RBF ')
plt.plot(X, y_lin, c='r', label='Modelo Lineal')
plt.plot(X, y_poly, c='b', label='Modelo Polinomial')
plt.xlabel('datos')
plt.ylabel('target')
plt.title('M. De Soport de Vectores')
plt.legend()
plt.show()

La gráfica que se muestra es la siguiente:

SVM_Regresión

En las primeras ocasiones que uno usa esta familia de algoritmo la confusión se genera con el concepto de “kernel”,  de modo no formal y burdo la idea de lo que hace el kernel es modificar el modo en el cual se logran definir los planos separadores o las fronteras entre clases de datos. Esto en palabras puede no parecer claro, pero el ejemplo natural para ejemplificar esto es generar una muestra de datos altamente no lineales y aplicar SVM con diferentes kernel lo cual mostrará como se comporta el algoritmo al cambiar el tipo de kernel. En el ejemplo hago una comparación con el kernel rbf, lineal, polinomial de grado 4 y 5, y el kernel sigmoid.

El código es el mismo para las cuatro estimaciones solo se cambiar una línea de código la cual para obtener las gráficas se comentan las no requeridas. El código es prácticamente una ligera modificación al ejemplo que se puede encontrar en scikit-learn.

El código para realizar este ejemplo es el siguiente:

#Código
#Module
print(__doc__)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVC
import time

start_time=time.time()
xx, yy = np.meshgrid(np.linspace(-3, 3, 1000),
 np.linspace(-3, 3, 1000))
np.random.seed(0)
X = np.random.randn(1000, 2)
Y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)


# fit the model

#clf =SVC(kernel='poly',degree=5)
clf =SVC(kernel='poly',degree=4)
#clf =SVC(kernel='linear')
#clf=SVC()
clf.fit(X, Y)

# plot the decision function for each datapoint on the grid
Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

plt.imshow(Z, interpolation='nearest',
 extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',
 origin='lower', cmap=plt.cm.PuOr_r)
contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,
 linetypes='--')
plt.scatter(X[:, 0], X[:, 1], s=30, c=Y, cmap=plt.cm.Paired)
plt.xticks(())
plt.yticks(())
plt.axis([-3, 3, -3, 3])
plt.title('Kernel Polinomio grado 4')
plt.show()

print(" The time required is %f seconds "%(time.time()-start_time))

Las gráficas que se obtienen son las siguientes:

SVM_Kernel_Lineal SVM_Kernel_Polinomial_degree4
SVM_Kernel_RBF

SVM_Kernel_Polinomial_degree5

Lo que se observa en las 4 gráficas es que el modelo con kernel rbf y un kernel polinomial de grado 4, parecen ser los dos mejores modelos para clasificar los datos en dos clases.Se observa que el modelo lineal no es para nada apropiado y que le cuesta trabajo al SVM con ese kernel poder separar los datos. El único modo que considero o creo que sirve para detectar si los datos son altamente lineales o altamente no lineales es hacer una exploración en los datos.

Un par de ejemplos.

Tomando datos de la plataforma http://datos.gob.mx/,  hago un ejemplo burdo,  ya que los datos que tomo son sobre los delitos por estado o Entidad Federativa en distintas clasificaciones de delitos, estos datos cubren desde el mes de Enero del 2001 hasta diciembre del 2006.  Lo que hago es medir la media de todos los delitos por cada uno de los registros de los datos, después teniendo la media defino una escala simple para asignar una variable indicadora a los estados que tienen mayor y menos tasa delictiva.

Con lo anterior defino un dos clases, debido a que no uso los datos para hacer un  análisis cuidadoso sino solo para comparar el uso del algoritmo, entonces omito la relevancia de tener las fechas de la medición por mes y año, tomo una muestra aleatoria de aproximadamente el 70% para entrenar el algoritmo y el resto para probar la eficiencia.

En resumen las etapas son las siguientes:

  1. Procesar y limpiar los datos.
  2. Asignar la variable de clasificación o indicadora.
  3. Entrenar el algoritmo y comparar con el conjunto de prueba.
  4. Comparar contra otro algoritmo, Naive Bayes.

 Primero recomiendo descargar los datos, ya que se cuenta con ellos mediar la media por cada fila y teniendo esa tabla la cargamos a python.

#Procesamiento
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
from pandas import Series

#Cargamos los datos
os.chdir('Directorio donde se pusieron los datos')

Datos=pd.read_csv('Archivos_con_los_datos.csv',sep=',',decimal=',')

Datos.head()
Datos.shape

#Exploración

Datos['Media'].hist()
plt.xlabel('Valores de la Media de la Delincuencia por Estado')
plt.ylabel('Valor de la Media')
plt.title('Histograma')
plt.show()

#Limpiamos los datos por los Missing values

sum(Datos.isnull().values.ravel())
#31

#Cleaning of data 
Datos1=Datos.dropna()
Datos1.shape 
#(2273,27)

S=Datos1['Media']
#Exploración 2
S2=[log(y) for y in S]

plt.hist(S2,col='green')
plt.xlabel('Valores de la Media de la Delincuencia por Estado')
plt.ylabel('Valor de la Media')
plt.title('Histograma')
plt.show()

#Missing Values
sum(Datos.isnull().values.ravel())
#31

#Cleaning of data 
Datos1=Datos.dropna()
Datos1.shape 
#(2273,27)

Series(S).mean()
#9.32


Hist_Delic_Mex

Histograma de las Medias de la Delincuencia por Estado

Hist_Del_Mex_Log

Histograma de los Logaritmos de las Medias de la Delincuencia por Estado

Después de explorar de modo breve como se comportan los datos ahora medimos a media creamos una nueva variable y le asignamos un valor dependiendo de la medición de la media. Agrego también los valores del logaritmo de la media ya que esta variable muestra que el comportamiento de las medias tienen una distribución de colas pesadas.

#Nueva variable
Datos1['Log_Media']=S2

Datos1.head()
Datos1.var()
 
#total data
 
len(Datos1[Datos1['Media']>9.3169]['Log_Media'])
#569
Datos1['Indicador']=0
Datos1.ix[Datos1.Media>9.3169,'Indicador']=1

Datos1.shape
#(2273,29)

Ahora preparo los datos para entrenar el algoritmo y aplico SVM, tanto con un kernel lineal como con uno no lineal.

# Preparación de datos y SVM
from sklearn.svm import SVC
from sklearn import svm
from sklearn.naive_bayes import GaussianNB
import random 

rows=random.sample(Datos1.index,1536)

Datos_train=Datos1.ix[rows]
Datos_test=Datos1.drop(rows)

type(Datos_train)
X=np.asarray(Datos_train.icol(range(4,26)))
Y=np.asarray(Datos_train.icol(28))
X_test=np.asarray(Datos_test.icol(range(4,26)))
Y_test=np.asarray(Datos_test.icol(28))
X_test

#Aplico los algoritmos de SVM

#Kernel no lineal
clf=SVC()
clf.fit(X,Y)
len(Y_test)
#737
P=clf.predict(X_test)
sum(P==Y_test)/737.0
#0.778

#Kernel Lineal
lin_clf=svm.LinearSVC()
lin_clf.fit(X,Y)
Pred=lin_clf.predict(X_test)
sum(Pred==Y_test)/737.0
#0.9701

#Naive Bayes
gnb = GaussianNB()
Prediction=gnb.fit(X,Y).predict(X_test)
(Prediction!=Y_test).sum()
(Prediction==Y_test).sum()
683.0/737.0
#0.926

Se observa que el porcentaje de eficiencia de los algoritmos es el siguiente, para SVM con un kernel no lineal es del 77.8%, para SVM con un kernel lineal es del 97.01% y para el algoritmo Naive Bayes es de 92.6%.

Entonces en resumen, para esta muestra de datos se obtienen que SVM con un kernel lineal es el más eficiente como método de clasificación.

La moraleja de siempre es usar más de un algoritmo para una muestra de datos, ya que ninguno es el “mejor” solamente es más adecuada para cada cierto tipo de datos.

Referencias:

1.- An Introduction to Statistical Learning

2.-The Elements of Statistical Learning

3.-Machine Learning: A Probabilistic Perspective

4.-Pattern Recognition  and Machine Learning

Encontrando grupos o análisis de Cluster

“The validation of clustering structures is the most difficult and frustrating part of cluster analysis. Without a strong effort in this direction, cluster analysis will remain a black art accessible only to those true believers who have experience and great courage”.. Jain and Dubes

La intención de esta entrada es complementar parte de los explicado en la entrada  “Análisis de Cluster, un ejemplo sencillo”. El código en esa entrada fue hecho en R project y expliqué pocos aspectos teóricos del tema.

En esta entrada intento cubrir los básico y dar más ejemplos de como usar esta técnica. El código  que comparto está hecho en Python  y poco a poco indico dónde y cómo descargar los datos de cada ejemplo.

La entrada anterior de esta categoría “Algo de sobre Python, análisis de datos y machine learning”, expliqué o di referencia sobre las librerías que usaría en el resto de entradas. Se debe de contar con ellas sino el código puede funcionar.

El ejemplo clásico “dato de Iris”.

Una muestra de datos clásica para analizar es sobre la morfología de las Iris setosa,  versicolor y virgínica, estos fueron inicialmente estudiados por R. Fisher aproximadamente en 1930.

Para ejemplificar la técnica uso dos librerías, mlpy y scikit-learn, en la última se cuenta con los datos de las Iris, así que solo se cargarán y se analizan primero gráficamente y luego se hace uso de un algoritmo estándar de análisis de Clusters.

Comentario: La ventaja de la librería Scikit-Learn es que pueden extraer datos de varios repositorios y además cuenta con ciertas bases de datos muestra, eso ayuda a practicar los algoritmos de ML.

#Cargamos los datos de Iris desde Scikit-learn

#Graficamos 

#Importamos las librerias necesarias

from matplotlib import pyplot as plt
from sklearn.datasets import load_iris

#Cargamos los datos y graficamos

datos=load_iris()
caract=datos.data
caract_names=datos.feature_names
tar=datos.target

# Graficamos los datos con colores distintos y tipo de marcas distintos

for t, marca,c in zip(xrange(3),">ox","rgb"):
 plt.scatter(caract[tar==t,0],caract[tar==t,1],marker=marca,c=c)
 
plt.show() 

Para correr el código de manera rápida se puede copiar en el IDLE que se descarga al instalar Python  y usar F5 o poner el código en un archivo “plano”, con terminación “.py” y correrlo en el cmd o consola. Los detalles se puede revisar en el libro digital “Learn Python the hard way”.

La mejor opción es hacer una carpeta y guardar los archivos .py , iniciar una sesión con IPython y ejecutar el código.

El programa anterior carga los datos y el gráfico. Siempre es recomendable al hacer análisis de datos explorar la información y la vía en general es revisar propiedades estadística, rango de datos, propiedades y por medio de algún gráfico visualizar la información.

La gráfica que obtenemos es la siguiente:

Cluster_1_1

Desde una postura un poco matemática se aprecia que los datos son puntos en el plano, en este caso en el primer cuadrante, al pintarlos de distinto color podemos apreciar su distribución y se aprecia como los datos están “mezclados”.

Ahora probamos correr un código donde exploraremos los datos mediante el algoritmo K-Medias, explico adelante en qué consiste este algoritmo.

#Importamos las librerias necesarias

from matplotlib import pyplot as plt
from sklearn.datasets import load_iris
import mlpy 
from sklearn.cluster import KMeans
 
#Cargamos los datos y graficamos

datos=load_iris()
dat=datos.data
caract_names=datos.feature_names
tar=datos.target

#Calculamos los cluster

cls, means, steps = mlpy.kmeans(dat, k=3, plus=True)

#steps
#Esta variable permite conocer los pasos que realizó el algoritmo para terminar

#Construimos las gráficas correspondiente

plt.subplot(2,1,1)
fig = plt.figure(1)
fig.suptitle("Ejemplo de k-medias",fontsize=15)
plot1 = plt.scatter(dat[:,0], dat[:,1], c=cls, alpha=0.75)
#Agregamos las Medias a las gráficas

plot2 = plt.scatter(means[:,0], means[:,1],c=[1,2,3], s=128, marker='d') 
#plt.show()

#Calculamos lo mismo mediante la librería scikit-lean

KM=KMeans(init='random',n_clusters=5).fit(dat)

#Extraemos las medias

L=KM.cluster_centers_

#Extraemos los valores usados para los calculos

Lab=KM.labels_

#Generamos las gráfica

plt.subplot(2,1,2)
fig1= plt.figure(1)
fig.suptitle("Ejemplo de k-medias",fontsize=15)
plot3= plt.scatter(dat[:,0], dat[:,1], c=Lab, alpha=0.75)

#Agregamos las Medias a las gráficas

plot4= plt.scatter(L[:,0], L[:,1],c=[1,2,3,4,5], s=128, marker='d')

#Mostramos la gráfica con los dos calculos

plt.show() 

El programa regresar una gráfica como la siguiente:

Cluster_4

 

Se aprecia en las gráficas que son los mismos datos pero pintados con distintas cantidad de colores, más aún en cada gráfica se encuentran 3 y 5 puntos, respectivamente, de tamaño mayor al resto.

Revisando el código anterior con cuidado, se observa que se hace uso de dos librerías mlpy y scikit-learn, en las cuales aplicamos el mismo algoritmo.

En la primera de las gráficas se aglomeran lo datos en 3 grupos y en la segunda en 5, que es lo que corresponde con los parámetros pedidos en las funciones mlpy.kmeans y KMeans respectivamente.

Entonces lo que se aprecia es que se estiman 3 y 5 clusters y además se inicializa el algoritmo de modo distinto (kmean++ y random) y para ver los detalles de las funciones se puede consultar en las páginas de referencia[1,2].

La inicialización del algoritmo es fundamental, ya que para procesar muchos datos la complejidad computacional se ve reflejada tanto de usos de memoria como de procesamientos, entonces es importante el método que se usa.

¿Qué con este ejemplo? Bueno, lo que se hace es explorar que se tienen 3 grupos de datos correspondientes a distintas Iris, luego se prueba explorar si el modelo de clusters describe los datos. Pero faltan muchos detalles, el primero que inmediatamente surge es confirmar cual es mejor modelo, el de 3 o 5 clusters, ¿cómo validar esto?, ¿qué usar para saber si el modelo es bueno?

Trato más adelante de contestar estas preguntas en otros ejemplos.

Pero lo anterior representa la idea escencial de análisis de Clusters, que es agrupar los datos en grupos similares[3].

El comentario anterior prácticamente dice todo, pero agrupar implica encontrar un algoritmo para hacerlo y por grupos similares implica saber cuando decir que dos cosas son similares, los dos detalles no son triviales.

Un poco de teoría

Por supuesto que lo que menciono como teoría es breve y para nada exhaustivo. Dejo algunas referencias que considero suficientes para revisar los aspectos teóricos, depende del gusto y formación [3,4,5].

Lo básico para entender los modelos de Cluster es tener la noción de métrica o distancia, con esa idea se define lo que es la matriz de similaridad y disimilaridad.

Métrica. En matemáticas esta noción permite el estudio de los espacios métricos, pero un ejemplo sencillo de métrica es el valor absoluto en los números reales, estos son los números que desde la primaria nos enseñan y los usamos para contar cosas con cierta precisión como 1.5 kilogramos de arroz, 2.15 metros, 355 ml, etc.

Observamos que si tenemos los número reales, el valor absoluto siempre es positivo o cero. Ejemplo, abs(2.15)= 2.15 y abs(-2.15)=2.15 y  para abs(0)=0.

Entonces podemos decir que si x es un número real, entonces abs(x) siempre es mayor o igual a cero.

Otra propiedad que tiene el valor absoluto es que el único número que tienen valor absoluto cero es el mismo cero, es decir; si abs(x)=0 entonces x=0.

La última implicación, para la gente que tiene poca familiaridad con matemáticas no resulta tan natural. Pero otra idea es pesar en caminatas, si del punto x al punto y se miden cierta cantidad de pasos y nos encontramos en el punto x. La distancia del punto x al mismo punto x debe ser cero y eso nos dice que sabiendo que nuestra posición inicial era el punto x, entonces estamos en el punto x. Medio en formalidad es, si distancia(x,y)=0, entonces x=y.

La siguiente idea se base es ver un triángulo rectángulo y medir sus lados con una regla, esto nos indicaría que la suma de los valores absolutos de su base y su altura es mayor o igual al valor absoluto de la hipotenusa. Si uno observa el teorema de Pitágoras, aprecia que a^2+b^2=c^2, pero eso no dice que la igual persista con los valores absolutos de a, b y c.

Entonces lo que se tiene es que abs(hipotenusa) es menor o igual a abs(altura)+abs(base). Esta propiedad se llama desigualdad del triángulo.

Lo que se tiene es que el valor abs() cumple con tres propiedades, que pueden ser vistas como la abstracción de medir con una regla líneas sobre una hoja en una mesa plana. Estas tres propiedades son lo que definen una métrica o distancia.

Recalco lo mencionado de “una mesa plana” por que si uno imagina medir rectar o segmentos de rectas sobre la superficie de un balón, la situación cambia drásticamente y da origen a cosas raras como triángulos que tienen como suma de sus ángulos más de 180°, eso es en general lo que se llama geometría no euclidiana o que el espacio no es euclidiano.

¿Por qué menciono eso de lo no-euclideano?

Ese punto en los algoritmos de clusters es importante, teóricamente para evitar complicaciones supone que los datos están en un espacio euclideano, pero en general muchos datos tienen comportamientos que implican ser modelados como no-euclideanos. Esto hace que se mejoren las técnicas y se estudien aspectos teóricos de los algoritmos[6].

Teniendo  en cuanta que la mayoría de datos que se analizan provienen de muchas variables, lo que tenemos es que las métricas o distancias que podemos usar con nuestros datos dependen de la naturaleza misma de los datos.

Es decir, si tenemos datos sobre 50 indicadores financieros medidos en un día, lo que podríamos considerar es usar la distancia Euclineana que es un modo de medir la distancia entre puntos de un plano o espacio, ya que posiblemente los valores de los datos se pueden modelar con lo números reales.

Pero no ocurre lo mismo si tenemos datos económicos y sociológicos de 100 países, ya que posiblemente varias variables serán categóricas u ordinales, para lo cual se debe de usar otro tipo de métricas.

 Similaridad y disimilaridad. La idea es tener una métrica y esto mide la cercanía o la distancia entre dos de nuestras variables, pero en la prática es difícil tener una distancia en el estricto sentido, sobre todo que satisfaga la desigualdad de triángulo. Lo que se hace es considerar que la métrica solo cumple las dos primeras propiedades, es decir, d(U,V)>0 y d(U,U)=0, donde U y V son variables.

Un poco de álgebra. En todos los algoritmos lo que se construye con las variables es una matriz de disimilaridad que tienen ciertas propiedades, se construye considerando la entrada Aij=d(Ui,Uj).

Los algoritmos piden que A sea una matriz simétrica o en caso de que A sea una matriz de similaridad se convierte en una de disimilaridad transformando la matriz por medio de una  función monótona no decreciente.

Si no se conoce mucho de matrices y sus propiedades, en una idea burda se pueden pensar en ellas como una tabla con números y su relevancia en los algoritmos  es que se sabe mucho de ellas para trabajar con operaciones matemáticas y representarlas en una computadora no es complicado. En especial en Python se cuenta con la librería Numpy que permite trabajar con matrices y en R project existe el tipo de dato matrix() para su manejo.

Métricas y tipo de variables. En resumen se pueden clasificar el tipo de variables en tres familias: cuantitativas, ordinales y categóricas [4].

Para cada tipo de variable se tienen varios tipo de métricas, entre las usuales son:

Para cuantitativas: Euclídea, Manhattan, L1, correlación estandarizada.

Para ordinales: Se convierten en cuantitativas, mediante un ajuste, si tenemos M valores posibles, de cambia cada entrada por (i-0.5)/M.

Para categóricas: la Discreta es usual.

Pero existen muchas métricas que pueden ser empleadas, dependiente del tipo de datos y análisis que se pretenda hacer, basta echar un vistazo a la variedad que existen.

Más aún, el concepto de métricas no restringe hacer uso de alguna métrica conocida, si uno define una regla para medir y se prueba que satisface por lo menos lo mínimo para obtener una medida de disimilaridad puede hacer uso de los algoritmos de análisis de Clusters.

Por su puesto que no es necesario aprenderse todas las métricas que se conocen, lo recomendable es explorar varias posibilidades de las métricas estándar por tipo de variable y probar con otras para explorar si existe algún beneficio en el análisis.

Sobre los Algoritmos. Los dos algoritmos más usados son agrupamiento jerárquico (Hierarchical Clustering) y K-Medias (K-means). Existen varias versiones y adecuaciones de ellos, sobre todo para casos cuando las condiciones implican que el modelo se encuentra en espacio no-euclideanos, donde los datos requieren cierto tipo de medida de disimilaridad o cuando se modelan datos de larga-escala[6].

Los detalles teóricos respecto a tiempo de convergencia y la descripción del algoritmo pueden ser consultadas en las referencias [4,5,6].

 Para ver como se comportan los algoritmos de un modo gráfico y con detalle de lo que pasa en cada etapa recomiendo consultar las notas de Andrew Moore[7].

Tomo dos imágenes sencillas de los algoritmos debidas a Toby Segaran[8].

HC

KM

El primero es una imagen para explicar como funciona el algoritmo jerárquico, la idea es que inicia con A y B, ve su disimilaridad y avanza el algoritmo para detectar que puede ser también agrupado C. Avanza un paso más y detecta que D y E estás “lejos”, por lo cual avanza y agrupa a estos últimos. Lo que se oculta la idea de un árbol en los datos, donde las primeras hojas con formadas por A y B, su padre es C y de manera independiente se encuentran las hojas D y E. Se ve en la siguiente figura:

AHC

El segundo algoritmo, K-media; en la imagen se inicia con dos puntos “aleatorios” en la primera etapa los cuales serían los centro de cada cluster. En la segunda etapa se calcula el centroide de los punto cercanos y se cambia el centro de los clusters al centroide calculado, se continua de este modo hasta que ya no hay variaciones en la posición de los clusters.

Los dos algoritmos tienen variaciones y definiciones de parámetros tratando de adecuarlos u optimizarlos ante cierto tipo de datos, en general el que k-medias tienen un tiempo menor de estimación con respecto al clusters jerárquicos.

Para mostrar lo simple que puede verse el algoritmo de k-medias describo los pasos desde el punto de vista de Machine Learning:

#Idea General del algoritmo
1)Inicialización
  -Se eligen la cantidad de clusters k
  -Se elige aleatóriamente k posiciones desde los datos de entrada
  -Se indica que el centro de los clusters tienen la posición de la media de los datos

2)Aprendizale
  -Se repiten los pasos:
      Para cada punto de los datos se hace:
       -Se calcula la distancia del punto al centro         del Cluster
       -Se aglomera los puntos con el centro del clu        ster más cercano
      Para cada Cluster
       -Se cambia la poción del cluster al centro al        centroide de los puntos agrupados o aglomera        dos
       -Se hacen cambios en el centro del Cluster ha        sta que ya no hay variación en el cambio de         posición
3)Tratamiento
      Para cada punto de prueba
       -Se calcula la distancia a cada cluster
       -Se asigna el punto al cluster con el cual ti        ene la distancia menor
      

 La visualización del algoritmo es tomado de la referencia [9], para ver en python como escribir estas indicaciones de manera básica solo se necesita la librería Numpy. El código es el siguiente:

#Implementación básica del Algoritmo en Python

from numpy import ones,argmin

#Calculo de distancia
dis=ones((1,self.nData))*sum((data-self.centre[0,:])**2,axis=1)

for i in range(self.k-1):
 dis=append(dis,ones((1,self.nData))*sum((data-self.centres[j-1,:])**2,axis=1),axis=0)

#Se identifican los clusters
cluster=dis.argmin(axis=0)
cluster=traspose(cluster*ones((1,self.nData)))

#Se actualizan los centros
for j in range(self.k):
 thisCluster=where(cluster==j,1,0)
 if sum(thisCluster)>0:
 self.centres[j,:]=sum(data*thisCluster,axis=0)/sum(thisCluster)

 La idea de compartir el anterior código y el esquema del algoritmo es mostrar como es sencillo el proceso del algoritmo e implementación, uno puede encontrar otras versiones en la referencia [8]. En los siguientes ejemplos no de hace uso del este código, se utiliza la librería scikit-learn, en particular el módulo llamado sklearn.cluster.

Clave: En programación no siempre es necesario inventar el hilo negro.

Ejemplos con datos

Los primeros dos ejemplos son variaciones a los ejemplos que se pueden encontrar en la página oficial de scikit-learn []. Lo que trato de hacer con estos dos ejemplos es explicar lo que hace el código y el algoritmo, pero sobre todo como son atractivos visualmente espero que esto motive un poco el interés sobre como usar estas técnicas.

En el segundo ejemplo es una variación y extensión un ejemplo de la referencia[10]. Para el cual explico como se procesan textos y cual es el objetivo al identificar clusters en este ejemplo, doy la referencia y el banco de bases de datos de donde puede descargar los archivos que se usan.

Ejemplo 1.Tomamos una imagen que vienen por default en scipy para hacer ejercicios de procesamiento de imágenes, es tan conocida que se llama Lena y corresponde a una fotografía de Lena Sjooblom una playmate de los años 70.

El código fue hecho por Vincent MichelAlexandre Gramfort, se puede consultar la versión original en example-cluster-plot-lean.

El código es el siguiente:

#Cargamos los datos de Iris desde Scikit-learn
#Manejo de imagenes y clusters

print(__doc__)
#Importamos las lib necesarias

import time as time
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from sklearn.feature_extraction.image import grid_to_graph
from sklearn.cluster import AgglomerativeClustering

####################################################

# Generando datos

lena = sp.misc.lena()

# Se reduce la resolucion de la imagen por un factor de 4

lena = lena[::2, ::2] + lena[1::2, ::2] + lena[::2, 1::2] + lena[1::2, 1::2]
X = np.reshape(lena, (-1, 1))

####################################################

# Se define la estructura de la conexion de la imagen, de los pixeles y sus vecinos

connectivity = grid_to_graph(*lena.shape)

####################################################

# Se calculo de Cluster
#Mandamos un mensaje de lo que está haciendo

print("Calculando la estructura de los clusters jerarquicos...\n")
print("No te desesperes ya viene la imagen")

st = time.time()
n_clusters1 = 20 # Numero de regiones estimadas en el cluster
ward1 = AgglomerativeClustering(n_clusters=n_clusters1,
 linkage='ward', connectivity=connectivity).fit(X)
label1 = np.reshape(ward1.labels_, lena.shape)
print("Tiempo transcurrido: ", time.time() - st)
print("Num de pixeles: ", label1.size)
print("Num de clusters: ", np.unique(label1).size)

st = time.time()
n_clusters2 = 15 # Numero de regiones estimadas
ward2 = AgglomerativeClustering(n_clusters=n_clusters2,
 linkage='ward', connectivity=connectivity).fit(X)
label2 = np.reshape(ward2.labels_, lena.shape)
print("Tiempo transcurrido: ", time.time() - st)
print("Num de pixeles: ", label2.size)
print("Num de clusters: ", np.unique(label2).size)

####################################################

# Se agregan los resultados a la imagen y se muestran las regiones estimadas


plt.figure('Imagen1',figsize=(5, 5))
plt.imshow(lena, cmap=plt.cm.gray)
for l in range(n_clusters1):
 plt.contour(label1 == l, contours=1,
 colors=[plt.cm.spectral(l / float(n_clusters1)), ])
plt.xticks(())
plt.yticks(())

plt.figure('Imagen2',figsize=(5, 5))
plt.imshow(lena, cmap=plt.cm.gray)
for l in range(n_clusters2):
 plt.contour(label2 == l, contours=1,
 colors=[plt.cm.spectral(l / float(n_clusters2)), ])
plt.xticks(())
plt.yticks(())

plt.show()

 Lo que obtenemos como resultado al correr el código son dos imágenes como las siguientes:

Lena1Lena2

Las imágenes muestra como se divide en 20 y 15 clusters, lo que uno debe de tener claro es que la imagen se traduce a información numérica y esta es la que se procesa al aplicar el algoritmo.

Ejemplo 2. Se hace una comparación con el modo en el cual se calculan los clusters de un conjunto de imágenes digitalizadas de lo números del 0 al 9, esta base de datos es muy conocida y se encuentras varias similares en otro repositorios. Los números se ven de esta forma:

14.4

El siguiente código hace uso de una herramienta de visualización. La idea es usar una “proyección”, con esto trato de decir que la información es originada de un espacio de muchas dimensiones y se visualiza en una imagen de 2 dimensiones, es decir; sobre un plano.

Algo así:

Proyecciones-conicas

Esta técnicas implica mencionar otros conceptos, pero trato de compartir este ejemplo por que permite ver como los clusters como una herramienta de clasificación y siempre “una imagen vale más de mil palabras”.

Este código fue hecho por Gael Varoquaux y se puede consultar en embebimientos digitales y clusters.

#Docstring en python
print(__doc__)
#Librerias usadas

from time import time
import numpy as np
from scipy import ndimage
from matplotlib import pyplot as plt
from sklearn import manifold, datasets
from sklearn.cluster import AgglomerativeClustering

#Se extraen lo datos
digits = datasets.load_digits(n_class=10)
X = digits.data
y = digits.target
n_samples, n_features = X.shape

np.random.seed(0)

def nudge_images(X, y):
 # Se varia la muestra de datos
 shift = lambda x: ndimage.shift(x.reshape((8, 8)),
 .3 * np.random.normal(size=2),
 mode='constant',
 ).ravel()
 X = np.concatenate([X, np.apply_along_axis(shift, 1, X)])
 Y = np.concatenate([y, y], axis=0)
 return X, Y


X, y = nudge_images(X, y)


#----------------------------------------------------------------------
# Visualizacion
#Se define esta funcion para simplificar el proceso de regreas 3 imagenes
def plot_clustering(X_red, X, labels, title=None):
 x_min, x_max = np.min(X_red, axis=0), np.max(X_red, axis=0)
 X_red = (X_red - x_min) / (x_max - x_min)

 plt.figure(figsize=(6, 4))
 for i in range(X_red.shape[0]):
 plt.text(X_red[i, 0], X_red[i, 1], str(y[i]),
 color=plt.cm.spectral(labels[i] / 10.),
 fontdict={'weight': 'bold', 'size': 9})

 plt.xticks([])
 plt.yticks([])
 if title is not None:
 plt.title(title, size=17)
 plt.axis('off')
 plt.tight_layout()

#----------------------------------------------------------------------
# Se hace un proceso geometrico llamado embebimiento para poder visualizar 
#el resultado del proceso de deteccion de clusters

print("Computo del embedimiento")
X_red = manifold.SpectralEmbedding(n_components=2).fit_transform(X)
print("Realizado.")


for linkage in ('ward', 'average', 'complete'):
 clustering = AgglomerativeClustering(linkage=linkage, n_clusters=10)
 t0 = time()
 clustering.fit(X_red)
 print("%s : %.2fs" % (linkage, time() - t0))

 plot_clustering(X_red, X, clustering.labels_, "%s linkage" % linkage)


plt.show()

 Al correr el código se obtienen 3 imágenes de clusters de los dígitos del 0 al 9, son hasta cierto punto parecidas pero si uno tienen calma y les dedica tiempo puedo observar ligeras diferencias.

Lo interesante en este ejemplo es ver como un conjunto de imágenes son clasificadas y agrupadas hasta cierto punto bien.

El método es cluster jerárquicos, pero una de las versiones de este algoritmo se llama “aglomerativo” que es como fue usado en la función AgglomerativeClustering. Lo que hace distinto entre cada uno de los algoritmo, average, single linkage y complete linkage es la definición de la medida de disimilaridad, se pueden ver detalles en la referencia [4].

Digitos1Digitos2Digitos3

Ejemplo 3.- Este último ejemplo requiere descargar una base de datos y revisar algunos aspectos del análisis de textos. Todo el proceso de realizará con scikit-learn y solo para ciertos detalles de hará uso de NLTK, que es el módulo por excelencia para hacer procesamiento de textos.

La idea de como procesar textos.

En general si uno escribe una oración, como: “Esta entrada ha sido muy larga, ¿por qué no logro terminarla?”.

Uno puede ver el texto como algo que hacemos de manera normal, pero lo que podemos pensar en cosas básicas como: cuántas palabras tiene, cuántos símbolos no son palabras, como medir que tan larga es una oración, etc.

Este tipo de cosas resultan ser calculables y visibles de manera rápida. Pero ,¿qué pasa con esto cuando tenemos toda una colección de textos?,¿cómo analizar sus relación entre palabras, frases  o textos completos?

Lo que a un matemático se le ocurriría inmediatamente es “medir la longitud de la frase” o definir un tipo de métrica. Ejemplo, “esntrada”, “entrada”, “etnrada” tienen longitud por letra 8, 7 y 7. Pero solo midiendo la longitud no podemos distinguir entre las últimas dos palabras.

Lo que se puede hacer es definir otra métrica como la ” métrica de Levenshtein“. El problema con solo medir este tipo de propiedades es que hasta cierto punto son locales, es decir; en cada palabra existe cierta medida y una frase u oración está formada por palabras. Eso no informa mucho respecto a como se comporta una frase completa.

Otra estrategia para analizar textos, no solo palabras; es construir lo que se llama ” bolsa de palabras”. Suponiendo que se tienen dos o tres oraciones  formadas por varias palabras, lo que se busca es medir la frecuencia de las palabras en cada frase. Ejemplo, las  “Esta entrada ha sido muy larga, ¿por qué no logro terminarla?” y “Esta entrada trata sobre técnicas de aprendizaje no supervizado”.

La bolsa de palabras estaría formada por todas la palabras diferentes de las dos frases y a cada frase se le asocia un vector de formado con ceros y unos que indican si apareció la palabra en la frase o no. Ejemplo para la primera oración forman una bolsa de palabras de 17 elementos, se tienen palabras en común, por lo cual el vector de la primera frase se vería como: (1,0,0,1,1,1,1,1,1,1,1,0,0,0,1,0,0).

Lo que observamos de esta representación es que conforme más frases con palabras diferentes la bolsa de palabras crece y la dimensión del vector también, con dimensión en este caso corresponde al número de entradas. Formalmente corresponde con la dimensión del espacio vectorial.

Otro beneficio de la representación por vectores es la posibilidad de definir una métrica mucho más geométrica, como la distancia euclideana.

Esto no es suficiente para preprocesar textos, otro elemento que se debe de considerar es eliminar las puntuaciones, ya que estos siguen siendo caracteres de la frases.

Un aspecto sutil, es que al establecer una métrica entre frases, puede resultar que una frase es repetida más de dos ocasiones, ejemplo “Cómo estas, cómo estas, cómo estas”. Dentro de la bolsa de palabras lo que uno obtendrá será un vector como el de la frase “Cómo estas” pero multiplicado por 3. Es decir, un vector de la forma (3,0,0,3,0,0,0,0,), entonces lo que se hace para evitar este tipo de anomalías es algo usual en álgebra lineal, normalizar los vectores. Es decir, se cambia el vector origina (3,0,0,3,0,0,0,0) dividido entre raíz de 2, lo cual hace que el vector sea igual en la norma al obtenido de (1,0,0,1,0,0,0,0). Así esas dos frases ya pueden considerarse iguales en el espacio de frases.

Lo último que se hace al analizar textos es eliminar las palabras más usuales, ya que estas aparecerán con mayor frecuencia en todas la frases y no permiten distinguir entre lo que realmente hace diferente o similar a una frase de otra. Además de que reduce la bolsa de palabras. Ya una vez que se “limpio” la frase resta analizar la relación de sus palabras o la frecuencia de las palabras respecto a la colección de textos, esto da una media de que tan relevante es la palabra. A este valor de la palabra se le llama Tf-idf.

Resumen:

1.-Se limpia la frase o texto de puntuaciones.

2.-Se elimina las palabras frecuentes del lenguaje

3.-Se construye el vector y se normalizan

4.-Se calcula el Tf-idf para las palabras que construyeron el espacio.

Todo lo anterior si bien mal explicado o incompleto es para tratar de dar una idea de lo que se hace con los textos. Al tener como medida de similaridad la distancia Euclideana podemos hacer uso del algoritmo de K-medias y elegir una cantidad de clusters para el algoritmo.

Para el ejemplo se hace uso de una base de datos llamada : 20newsgroup. Se puede descargar de varios modos, directamente desde el servidor http://qwone.com/~jason/20Newsgroups/. Otra opción es descargarlo desde un repositorio que cuenta con muchas bases de datos y se encuentra en http://mlcomp.org/datasets/379. La opción que considero más simple, pero se requiere un poco de paciencia es descargarla haciendo uso del módulo scikit-learn, mediante la función sklearn.datasets, las instrucciones se pueden leer aquí: http://scikit-learn.org/dev/datasets/index.html#datasets

El código se puede descargar desde https://github.com/luispedro/BuildingMachineLearningSystemsWithPython, el código se encuentra en el capítulo 3.

Las instrucciones para replicar el ejemplo:

1.-Descargar la base de datos

#Correr este código en Python o en IPython

from sklearn.datasets import fetch_20newsgroups
>>> newsgroups_train = fetch_20newsgroups(subset='train')

>>> from pprint import pprint
>>> pprint(list(newsgroups_train.target_names))
['alt.atheism',
 'comp.graphics',
 'comp.os.ms-windows.misc',
 'comp.sys.ibm.pc.hardware',
 'comp.sys.mac.hardware',
 'comp.windows.x',
 'misc.forsale',
 'rec.autos',
 'rec.motorcycles',
 'rec.sport.baseball',
 'rec.sport.hockey',
 'sci.crypt',
 'sci.electronics',
 'sci.med',
 'sci.space',
 'soc.religion.christian',
 'talk.politics.guns',
 'talk.politics.mideast',
 'talk.politics.misc',
 'talk.religion.misc']

2.-Guardar el código y correrlos desde:

-La consola estando en el directorio donde esta el código poniendo 
C:\......\>python nombredelarchivo.py
-Iniciar sesión en IPython y correr el código
In [1]%run python nobredelarchivo.py

El código del ejemplo es el siguiente, dejo los datos de origen de los creadores del código y solo agregué comentarios y explicaciones de lo que hace.

#Dejo los comentarios de origen respetando el trabajo de los autores


# This code is supporting material for the book
# Building Machine Learning Systems with Python
# by Willi Richert and Luis Pedro Coelho
# published by PACKT Publishing
#
# It is made available under the MIT License


#Cargamos los módulos requeridos

import sklearn.datasets
import scipy as sp


#Este será el post de prueba

new_post = \
 """Disk drive problems. Hi, I have a problem with my hard disk.
After 1 year it is working only sporadically now.
I tried to format it, but now it doesn't boot any more.
Any ideas? Thanks.
"""
print '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n'
print "Este sera nuestro poste de prueba\n"

print("""\
Dear reader of the 1st edition of 'Building Machine Learning Systems with Python'!
For the 2nd edition we introduced a couple of changes that will result into
results that differ from the results in the 1st edition.
E.g. we now fully rely on scikit's fetch_20newsgroups() instead of requiring
you to download the data manually from MLCOMP.
If you have any questions, please ask at http://www.twotoreal.com
""")

print ('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
#Cargamos la base de datos

all_data = sklearn.datasets.fetch_20newsgroups(subset="all")

print("Numero total de posts: %i" % len(all_data.filenames))
# Number of total posts: 18846

#Los grupos que se analizarán

groups = [
 'comp.graphics', 'comp.os.ms-windows.misc', 'comp.sys.ibm.pc.hardware',
 'comp.sys.mac.hardware', 'comp.windows.x', 'sci.space']
#Se elige la parte de la base a la cual se aplicará el algoritmo

train_data = sklearn.datasets.fetch_20newsgroups(subset="train",
 categories=groups)

print ("Numero de post de entrenamiento en los grupos de tecnologia :", len(train_data.filenames))

# Number of training posts in tech groups: 3529

print ("Se tomaran 60 Clusters!!")
labels = train_data.target

num_clusters = 60 # sp.unique(labels).shape[0]

#Se inicia el procesamiento de los datos

import nltk.stem
english_stemmer = nltk.stem.SnowballStemmer('english')

from sklearn.feature_extraction.text import TfidfVectorizer

#Se define una clase para asignar el Tfdf

class StemmedTfidfVectorizer(TfidfVectorizer):

 def build_analyzer(self):
 analyzer = super(TfidfVectorizer, self).build_analyzer()
 return lambda doc: (english_stemmer.stem(w) for w in analyzer(doc))

#Se continua con el preprocesamiento de los datos

vectorizer = StemmedTfidfVectorizer(min_df=10, max_df=0.5,
 stop_words='english', decode_error='ignore'
 )

#Se contruyen el espacio de vectores que representan a los post que se analizan

vectorized = vectorizer.fit_transform(train_data.data)
num_samples, num_features = vectorized.shape

#Se presenta la cantidad de muestra

print("#Muestra: %d, #Caracteristicas: %d" % (num_samples, num_features))

# samples: 3529, #features: 4712

#Se importa la libreria para aplicar k-medias

from sklearn.cluster import KMeans

km = KMeans(n_clusters=num_clusters, n_init=1, verbose=1, random_state=3)
clustered = km.fit(vectorized)


print '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'

print("km.labels_=%s" % km.labels_)

# km.labels_=[ 6 34 22 ..., 2 21 26]


print("km.labels_.shape=%s" % km.labels_.shape)
# km.labels_.shape=3529

#Se cargan el submódulo de métricas para hacer la medición de similaridad

from sklearn import metrics

print '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'

print 'Metricas'

print("Homogeneity: %0.3f" % metrics.homogeneity_score(labels, km.labels_))
# Homogeneity: 0.400
print("Completeness: %0.3f" % metrics.completeness_score(labels, km.labels_))
# Completeness: 0.206
print("V-measure: %0.3f" % metrics.v_measure_score(labels, km.labels_))
# V-measure: 0.272
print("Adjusted Rand Index: %0.3f" %
 metrics.adjusted_rand_score(labels, km.labels_))
# Adjusted Rand Index: 0.064
print("Adjusted Mutual Information: %0.3f" %
 metrics.adjusted_mutual_info_score(labels, km.labels_))
# Adjusted Mutual Information: 0.197
print(("Silhouette Coefficient: %0.3f" %
 metrics.silhouette_score(vectorized, labels, sample_size=1000)))
# Silhouette Coefficient: 0.006

#Se procesa el post de prueba

new_post_vec = vectorizer.transform([new_post])
new_post_label = km.predict(new_post_vec)[0]

similar_indices = (km.labels_ == new_post_label).nonzero()[0]

similar = []
for i in similar_indices:
 dist = sp.linalg.norm((new_post_vec - vectorized[i]).toarray())
 similar.append((dist, train_data.data[i]))

similar = sorted(similar)
print '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
print("Count similar: %i" % len(similar))

print '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
#Se eligen los tres más similares
show_at_1 = similar[0]
show_at_2 = similar[int(len(similar) / 10)]
show_at_3 = similar[int(len(similar) / 2)]

#Se imprimen

print("=========================== #1 =========================")
print(show_at_1)
print()

print("=========================== #2 =========================")
print(show_at_2)
print()

print("=========================== #3 =========================")
print(show_at_3)

El código anterior parece raro o feo, a primera vista pero lo que hace es bastante interesante.

En primera se toman 3529 elementos de muestra para analizar, de entre ellos se construye un espacio con los vectores asignados a cada post, claro que para llegar a eso primero se procesan los textos.

Teniendo el espacio vectorial se procede a medir la distancia, esto permite la construcción de una matriz de términos, que permite aplicar el método de k-medias sobre los post vectorizados.

Teniendo esto se hacen mediciones de métricas usuales sobre los clusters.  Sobre estos indices o medidas, no menciono nada, pero se pueden consultar en la referencia [3].

Cuando los datos se han dividido por clusters, entonces se procede a calcular o estimar en cual es el cluster que le corresponde al post y se eligen los 3 post que tienen menor similaridad o distancia euclideana.Termina el programa presentando los 3 post más cercanos.

Espero este último ejemplo no desanime a la primera lectura, es menos visual que los otros dos ejemplos previos pero que es bastante interesante, no tan solo por que procesa una base considerablemente grande, sino por que se muestra como emplear el algoritmo de clusters en una actividad que no suele mencionarse en los textos con perfil muy estadístico.

Este tipo de procedimiento, es decir; procesar los datos y luego asignar un tipo de medida, puede ser llevado a la práctica con otro tipo de información que no solo sea texto. El más interesante es quizás el analizar base de datos de imágenes y comparar los resultados con otras técnicas de clasificación.

Referencias.

1.-http://en.wikipedia.org/wiki/K-means%2B%2B

2.-http://www.autonlab.org/tutorials/kmeans11.pdf

3.-Capítulo 25, http://www.cs.ubc.ca/~murphyk/MLbook/

4.-Capítulo 13 y 14, http://statweb.stanford.edu/~tibs/ElemStatLearn/

5.-Capítulo 9, http://research.microsoft.com/en-us/um/people/cmbishop/prml/index.htm

6.- Capitulo 7, http://www.mmds.org/

7.-http://www.autonlab.org/tutorials/kmeans11.pdf

8.-http://shop.oreilly.com/product/9780596529321.do

9.-Capitulo9,http://www.amazon.com/Machine-Learning-Algorithmic-Perspective-Recognition/dp/1420067184

10.Capitulo3y4, http://shop.oreilly.com/product/9781782161400.do

11.-Página Oficial de Scikit-Learn

Sistemas de recomendación

La entrada es breve y la divido en dos partes, un poco la teoría  y un poco la parte de ejemplos. Al final de la entrada dejo las referencias.

Cabe mencionar que en general ya no es solo un tópico el tema de los sistemas de recomendación, existen ahora cursos  relativamente amplios respecto al tema, el recomendado es el curso de la universidad de Minessota y que cuentan con el laboratorio de computo GroupLens en el cual se estudian técnicas de sistemas de recomendación.

El impacto en general para todos es conocido, uno puede ver que youtube recomienda lista de canciones, que Amazon recomienda libros, que Coursera recomienda cursos y el ahora desaparecido Grooveshark recomendaba géneros y canciones, lo que ahora hace Spotify que entre sus títulos publicitarios dice “Déjanos entenderte”.

En general la frase de Spotify es lo que significan los sistemas de recomendación, tratar de entender las preferencias de los usuarios para recomendarte algún producto o servicio.

El caso más conocido es el trabajo realizado por dos equipos en el premio de Netflix, los equipos tomaron el nombre de Pragmatic Chaos y ayudaron a mejorar sus sistema de recomendación de la empresa, para más detalles pueden consultar las referencias [1,2,3,4].

Un breve recorrido por la teoría.

Esta información la tomo del curso “Mining of Massive datasets” de la universidad de Stanford, el cual es de libre descarga y la otra fuente recomendables, como ya mencioné, son las notas del curso de la Universidad de Minessota. Dejo esta dos fuentes en las referencias.

La clasificación de los sistemas de recomendación en general son dos:

  • Sistemas basados en contenido
  • Sistemas de Filtro colaborativo, en ingles es “Colaborative-Filtering”

El primero trata de recomendar tomando información del producto , sus propiedades y como este es “similar” a otro. El segundo tipo se basa en entender la relación entre el usuario y los productos.

En general el modelo para los dos tipos de sistemas es considerar una matriz de usuario y productos. Ejemplo, se tienen una lista de películas  y una lista de personas, no todas la películas fueron vistas por todos los usuarios pero más aún a cada usuario se le puede pedir su calificación de la película del 1 al 5. Así se puede ver la calificación de los usarios vs películas, donde las filas son los usuarios y las columnas las películas.

Lo delicado con esta matriz es el llenado, ya que en caso de tener muchos productos y una cantidad mucho menos de usuarios, la matriz estará casi vacía. Así que se debe de definir el tipo de estrategia para identificar o estimar las posibles preferencias de los usuarios con respecto a los productos. Esto da pie a la división en el tipo de sistemas de recomendación.

Por terminar….

1.-http://www.wired.com/2009/09/how-the-netflix-prize-was-won/

2.-http://www2.research.att.com/~volinsky/netflix/bpc.html

3.-http://www2.research.att.com/~volinsky/netflix/

4.-http://www.wired.com/2009/09/bellkors-pragmatic-chaos-wins-1-million-netflix-prize/

5.-https://www.coursera.org/learn/recommender-systems

6.-http://www.netflixprize.com/assets/ProgressPrize2008_BigChaos.pdf

7.-http://grouplens.org/

Algo sobre Python, análisis de datos y machine learning.

Sobre Python

Python es un lenguaje de programación interpretado multiparadigma, es decir; soporta hacer programación orientada a objetos y programación imperativa. Su creador fue Guido Van Rossum y debido a que es de código libre la comunidad de desarrolladores han creado librerías o módulos para hacer casi cualquier cosa.

En general la gente que hace análisis de datos o estadística conoce bien R project, SAS o SPSS. Pero pocos se han acercado a Python y sus librerías para análisis de datos, las ventajas pueden ser cuestionables con respecto a software que específicamente fueron diseñados con una perspectiva estadística, pero sin duda la potencia y calidad de librerías es muy buena. Principalmente para hacer uso de algoritmos de Machine Learning las librerías en Python son mejores que las de R project, sobre todo las técnicas de Deep Learning.

Quizás el mejor candidato para comparar el uso de Python en el análisis de datos es R project, por ser software libre y ser actualmente de alta demanda y atracción en la ciencia de datos. Sin embargo, Python siendo realmente un lenguaje de programación, no como R project, no se limita a ser usado  solamente para analizar datos, sino bien puede ser parte de un sistema o para desarrollar un proyecto completo. 

En caso de que resulte complicado hacer implementación de alguna técnica estadística en Python, se puede enviar un Script a R o manipular los datos desde R y extraerlos a Python por medio de la librería rpy2.

Yo siempre recomiendo para aprender a programar Python, por cosas muy sencillas, es fácil de aprender, fácil de leer y cuenta con librerías para hacer de todo. La curva de aprendizaje es muy corta y en muy poco tiempo uno puede estar programando cosas no tan triviales.

En caso de que se use Windows como sistema operativo, que la mayoría lo usa; se debe de instalar y revisar como hacer instalación de módulos, es bastante fácil y en general se hace desde la consola del sistema con el módulo pip. Si se usa Linux, ya no se tienen ningún problema, está por default, solo basta abrir la terminal del sistema para probarlo.

La mejor referencias para revisar información sobre python, tanto para instalar y conocer el lenguaje, es su página oficial:

https://www.python.org/

Para aprender desde cero a programar en este lenguaje existen muchas páginas y libros, pero creo que solo con ejemplos y ejercicios se aprende, recomiendo seguir y revisar el siguiente libro online:

http://learnpythonthehardway.org/book/

Comentario: cada ocasión que pongo un [] con un número me refiero al número de la lista de referencias.

Sobre Machine Learning y data analysis en Python

Considerando que se tiene instalado Python, las librerías necesarias para hacer los ejemplos que comentaré en las entradas de esta categoría son:

  1. Numpy
  2. SciPy
  3. Matplotlib
  4. Pandas
  5. mlpy
  6. scikit-learn
  7. Pybrain
  8. IPython
  9. NLTK

Por supuesto que hay más librerías necesarias, pero las comento conforme se haga uso de ellas.

Algunas de las convenciones que se hace al programar con las anteriores librerías, sobre todo con Numpy, SciPy , Pandas y Matplotlib, son las siguientes:

#Convenciones en los programas para hacer uso de ciertas librerías
import numpy as np
import pandas as pd
import scipy as sp
import matplotlib.pyplot as plt

 Lo único que dice el código anterior es que se hace uso de nombres cortos y que ahora son estándar en los programas que requieren estas librerías.

Suponiendo que ya se tiene instalado python y las librerías, uno puede revisar en la documentación del lenguaje en la página para ver como se hace uso con: números, variables boolenas, cadenas, listas, diccionarios, loops, tuplas, conjuntos, funciones y clases.

Principalmente el manejo de listas y diccionarios es de uso frecuente al analizar datos, ya que en general se usan para implementar algoritmos sobre un conjunto de datos, y cada uno de estas estructuras de datos tienen sus operaciones o métodos.

Si haz leído hasta este punto y te desanima un poco el aprender a instalar cosas, revisar tanta documentación para solo hacer uso de un software. Bueno, una opción es instalar un paquete con todas la librerías más usadas. Pensando en que mucha gente usa Windows y que la gente que usa Mac o Linux tienen familiaridad con instalación de software, entonces explico cuales son las opciones que conozco.

Las opciones son variadas y depende de las preferencias:

1.-Anaconda.

2.-WinonPyth

3.-Conopy

Cada una de estas instalaciones contiene todas la librerías que uso en los ejemplos, con excepción de pybrain y mlpy. Basta con seguir las instrucciones de instalación que son cosa de dar “click”.

Yo recomiendo Anaconda, es la que he estado usando recientemente y creo que está bien en general. Estos paquetes contienen más herramientas de las que menciono en la entrada, se puede revisar dependiendo del interés de cada uno para qué sirve el resto de librerías.

Un poco de NumPy, SciPy y Matplotlib

En el resto de la entrada trataré de explicar algunas operaciones y usos de las librerías Numpy, Scipy y Matplotlib.

Prácticamente son la base para muchas librerías, Numpy permite manipular arreglos o matrices de datos, si se tiene un poco de conocimiento de álgebra lineal resultará natural conocer las operaciones y funciones de dicha librería. Si no se cuenta con conocimiento de álgebra lineal, recomiendo revisar cualquier texto que prefieran ya que el perfil del texto va desde muy abstracto hasta muy aplicado y depende de la formación e interes de cada persona.

Pero es altamente recomendable revisar conceptos de estos temas lo mejor posible, ya que las matrices son el centro de muchos algoritmos de Machine Learning a todos los niveles.

SciPy, concentra varios procesos, algoritmos o técnicas de manipulación numérica; esto va desde el cálculo de la transformada rápida de Fourier, optimización, procesamiento de señales e imágenes.

Si nunca haz escuchado cosas de estos temas, no te desesperes ni mucho menos tienes que saber todo sobre eso, se aprende con la práctica y la formación, hay gente con carreras completas en ingeniería que no entendió la relevancia del análisis de Fourier  u optimización en sus cursos. En buena medida depende del problema que se aborde para hacer uso de alguna herramienta de Scipy.

Por último Matplotlib, es un módulo para gráficar, existe otro módulo para realizar gráficas que es actualmente muy usado, ggplot2. Este último es muy popular en la comunidad de usuarios de R project pero existe también la versión del módulo en Python y depende la preferencia en cuanto al tipo de gráficos al elegir entre Matplotlib o ggplot2.

El nivel de gráficos es bastante bueno y en ocasiones es tan importante hacer una buena gráfica ya que esto clarifica lo que se hace en el análisis de datos o como se presentan los resultados.

En los últimos años de han desarrollado muchos proyectos de visualización de datos y es prácticamente un campo de investigación, ya que al analizar muchos datos es fundamental tener un modo agradable de visualizar lo resultados de los análisis o investigaciones, de modo tal que sea interactivo e informativo [1].

¿Cómo correr el código en sus computadoras?

Los ejemplos los corrí en una sistema Windows con una versión del sistema 8.1 y con la versión de Python 2.7.9, lamento no comentar respecto a como correrlos en Linux, pero es mucho más simple ya que solo se debe de agregar la ruta del directorio en un script y debe de correr el código sin problema alguno.

Entonces pensando en que hay más usuario de Windows que de Linux, los detalles para correr el código que comentaré requiere probar que está bien instalado Python , configurar la variable PATH para correr Python desde la consola cmd o desde powershell y en caso de que se haga uso de Notepad++, se puede configurar para enviar el código a la consola y correrlos en el IDLE.

Creo que la mejor recomendación es empezar a usar Ipython y cualquier editor de texto, puede ser simplemente el notepad, pero se pueden revisar las recomendaciones en el libro online que mencioné anteriormente para ver como correr un scrip de python y como interactuar con el intérprete. Si presionan en la palabra “interprete”, podrán ver una versión online en la página oficial de python.

El párrafo anterior desalienta un poco cuando uno usa R project, por que solo basta instalar R, instalar las librerías y listo. Pero usar Python da oportunidad de conocer mejor como funciona aspecto del sistema, de la creación de una aplicación y de una que otra cosa que termina sirviendo posteriormente.

Lo que suele pasar con los usuarios de R es que pocas ocasiones desarrollan Script de R y los corren en la consola del sistema, el cmd en Windows. Esto cuando uno está iniciando a programar es raro, pero es fundamental para poder ponerse hacer cosas más interesantes.

Algunas aclaraciones, Notepad++ es un editor de texto más amigable y con más detalles y herramientas que el notepad normal, es de libre descarga y entre los lenguajes que soporta está Python. El IDLE que vienen en la descarga usual de Python es sencillo y lo único que se debe de revisar es como crear un archivo para ejecutar el código, que es cosa simple. Por supuesto que existen otros editores de texto y otros IDLE, depende del gusto de cada persona, pero la idea es la misma, de hecho directamente los IDLE no requieren otro editor como  Notepad++, en mi caso es por que me gusta más usar este último que escribir código directamente en el IDLE. Por último, lo recomendable es usar Ipython y una editor de texto eso basta para tener un entorno de trabajo con todos los recursos para trabajar con Python.

Un poco de Numpy

Esta librería permite la manipulación de arreglos o matrices de datos, en consecuencia también el manejo de vectores que son matrices de un solo renglón. Permite hacer operaciones aritméticas sobre todos los elementos de la matriz o elegir entre ellos bajo condiciones pedidas.

Haciendo un poco de memoria sobre las operaciones con matrices en álgebra lineal, en general uno estudia la dimensión de la matriz, la forma, el tipo de entradas que tiene, la multiplicación de matrices por un escalar, la suma y resta de matrices, el cálculo de la inversa de una matriz, el producto de matrices, el producto de una matriz y un vector, el producto interno de dos vectores o de las columnas de una matriz.

Este tipo de operaciones algebráicas se realizan con Numpy; algo extra es la manipulación de las entradas cuando se cuentan con datos NaN o eligiendo aquellas entradas que son mayores o menores algún valor, los que han usado Ovecta o Matlab.

En el código doy algunos ejemplos sencillos.

#Numpy
import numpy as np
#Creamos una matriz
a=np.array([0,1,2,3,4,5])
print 'Presenta la matriz'
print a
print 'Dimension de la matriz'
print a.ndim
#Vemos sus dimesión
print 'Forma de la matriz'
print a.shape
#Vemos su forma, el número de renglones y columnas
print 'Multiplica cada entrada por 2'
print a*2 
#Multiplicamos cada entrada por 2
print 'Eleva a la potencia 2 cada entrada'
print a**2
#Elevamos a la segunda potencia cada entrada
print 'Cambiando la forma del la matriz'
print a.reshape((3,2))
print 'Se imprime cada uno de sus elementos'
print a[0],a[1],a[2]
a[0]=5
print 'Se cambia la primera entrada a[0]=5'
print a
#Creamos una nueva matriz para de más dimensiones
b=np.array([[0,1,2],[3,4,5]])
print 'Se presenta b'
print b
print 'La forma de b'
print b.shape
print 'Se multiplica cada entrada por 3'
print b*3
print 'Se eleva cada entrada al cuadrado'
print b**2
#Se eliminan el vector y la matriz
del(a,b)


#Creacion de matrices especiales

a=np.zeros((3,4))
print 'Se presenta la matriz nula de 3X4'
print a
b=np.ones((2,2))

print 'Se presenta la matriz de unos'
print b
#Se creo la matriz con un valor igual para todas las entradas
c=np.full((3,3),3.5)
print 'Se presenta la matriz con valores iguales en todas las entradas'
print c

#Matriz identidad

d=np.eye(5)

print 'Se presenta la matriz identida de dimension 5'
print d

#Matriz con valores aleatorios

e=np.random.random((7,7))
print 'Se preseta la matri con valores aleatorios en cada entrada'
print e
#Eliminamos las matrices
del(a,b,c,d,e)

#Algunos detalles con los valores de las matrices y asignaciones

a=np.array([[1,2,3],[4,5,6],[7,8,9]])
print 'Se presenta la nueva matriz a'
print a
print 'Se presenta el valor de sus colnas'
print a[:,2]
#Se presentan las dos primeras filas
print 'Se presentan las dos primeras filas'
print a[:2,:]
#Creamos una nueva matriz con las primeras 4 entradas de a

b=a[:2,:2]
print 'Sub matriz de a'
print b
print 'modificamos el valor de b y modifica el de a'
b[0][0]=65
print b
print 'Modifica a'
print a
#esto tiene que ver con la asignación de memoria y de valores por python
del(a,b)
#Operaciones matematicas
a=np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=np.float64)
b=np.array([[1,0,1],[2,6,1],[1,1,1]], dtype=np.float64)
#Suma de matrices
print a+b
print np.add(a,b)
#Resta de matrices
print a-b
print np.subtract(a,b)

#Multiplicacion
print a*b
print np.multiply(a,b) 

#division
print b/a 
print np.divide(b,a)

#Raiz cuadrada

print np.sqrt(a) 

#Producto interno de vectores y de vectores con matrices

c=np.array([1,0,1])

print 'Produto interno por (1,0,1)'
print a.dot(c)
print 'Producto punto de matrices'
print a.dot(b)

print 'Valores de b'
print b
print 'elecci󮠤e valores mayores a 1'
print b[b>1]
print b>1
print 'Generamos otra matriz con un elemento NAN'
c=np.array([1,2,np.NAN,3,4])
print 'Presentamos su valor'
print c
print 'Identificamos los valores que son NAN'
print np.isnan(c)
print 'Elegimos los valores que no son NAN'
print c[~np.isnan(c)]

Puede haber evitado poner tanto comentario en pantalla o estar enviando “print” en cada orden y comentario, pero los dejo así sobre todo para gente que recién inicia a manejar Python y le pueda causar confusión solo ver números y no saber que cosa es lo que está haciendo.

Si se pone atención el objeto sobre el cual se trabaja siempre es np y de él se toman lo necesitado para crear el tipo de matrices que se desea.

Algunas cosas que no comenté son las operaciones como suma de valores en la columna o de valores en la fila, la transpuesta y el asignar a una matriz para que tengas valores enteros o flotantes. Esto es sencillo y se puede consultar en la documentación de NumPy.

Algo de SciPy y Matplotlib

Como antes mencioné SciPy es una caja de herramientas y Matplotlib es una librería para graficar. Dos ejemplo son los siguientes:

#Dos ejemplos de Scipy
#Cargamos los módulos necesarios
import numpy as np
from scipy import signal, misc
import matplotlib.pyplot as plt

#Se revisará una imagen como dato

image = misc.lena().astype(np.float32)
derfilt = np.array([1.0, -2, 1.0], dtype=np.float32)
ck = signal.cspline2d(image, 8.0)
deriv = (signal.sepfir2d(ck, derfilt, [1]) +
 signal.sepfir2d(ck, [1], derfilt))

#Se calcula el Lapaciano

laplacian = np.array([[0,1,0], [1,-4,1], [0,1,0]], dtype=np.float32)
deriv2 = signal.convolve2d(ck,laplacian,mode='same',boundary='symm')

#Se presenta la imagen
plt.figure()
plt.imshow(image)
plt.gray()
plt.title('Original image')
plt.show()

#Se presenta la imagen modificada por el filtro y la señal gaussiana

image = misc.lena()
w = signal.gaussian(50, 5.0)
image_new = signal.sepfir2d(image, w, w)
plt.figure()
plt.imshow(image_new)
plt.gray()
plt.title('Filtered image')
plt.show()

Las imágenes  que obtenemos es:

figura_1

 

Este es un ejemplo clásico de SciPy, la siguiente foto, generada por las últimas líneas del programa anterior no regresan la siguiente modificación a la imagen. Qué lo que se hace es extraer como señal gaussiana o filtrar la imagen original, esto lo pueden consultar en la página de SciPy o en textos de procesamiento de imágenes.

figure_2

 

Otro ejemplo estándar para usar SciPy es en el análisis de señales.

#Análisis de señales
from numpy import arange, sin, pi, random, array
x = arange(0, 6e-2, 6e-2 / 30)
A, k, theta = 10, 1.0 / 3e-2, pi / 6
y_true = A * sin(2 * pi * k * x + theta)
y_meas = y_true + 2*random.randn(len(x))

def residuals(p, y, x):
 A, k, theta = p
 err = y - A * sin(2 * pi * k * x + theta)
 return err

def peval(x, p):
 return p[0] * sin(2 * pi * p[1] * x + p[2])

p0 = [8, 1 / 2.3e-2, pi / 3]
print(array(p0))


from scipy.optimize import leastsq
plsq = leastsq(residuals, p0, args=(y_meas, x))
print(plsq[0])


print(array([A, k, theta]))


import matplotlib.pyplot as plt
plt.plot(x, peval(x, plsq[0]),x,y_meas,'o',x,y_true)
plt.title('Datos con Ruido')
plt.legend(['Estimacion', 'Ruido', 'True'])
plt.show()

figure_3

Observamos que todo el proceso anterior es el que se realiza para analizar una señal, se considera el ruido, se estima una curva y se analiza con respecto a los datos originales.

Los dos códigos anteriores son ejemplos de como se usa SciPy en dos contextos distintos, uno con números y otro con imágenes. Se puede consultar muchos ejemplos en la página oficial de SciPy y Matplotlib.

El último ejemplo es tomado del texto de “Building Machine Learning Systems with Python” y los datos se pueden descargar desde GitHub. La idea de lo siguiente es graficar el cruce entre datos de tráfico web y semanas, pero además agregar curvas ajustadas a los datos.

#Se simplifica el código presentado en el texto
import os
from utils import DATA_DIR
import scipy as sp
import matplotlib.pyplot as plt

 
sp.random.seed(3) 
data = sp.genfromtxt(os.path.join(DATA_DIR, "web_traffic.tsv"), delimiter="\t")
print(data[:10])
print(data.shape)

#Se separan los datos en x e y
x = data[:,0]
y = data[:,1]
print("Entradas no validas:", sp.sum(sp.isnan(y)))

#Se limpian los datos
x = x[~sp.isnan(y)]
y = y[~sp.isnan(y)]

#Modelo 1
fp1= sp.polyfit(x, y, 1)
f1 = sp.poly1d(fp1)
fx = sp.linspace(0,x[-1], 1000)
#Modelo 2
f2p = sp.polyfit(x, y, 2)
f2 = sp.poly1d(f2p)

#Gráfica de los modelos y los datos 
plt.scatter(x,y)
plt.title("Trafico Web en un Mes")
plt.xlabel("Tiempo")
plt.ylabel("Hits/hora")
plt.xticks([w*7*24 for w in range(10)],['Semana %i'%w for w in range(10)])
#Se agregan los ajustes
plt.plot(fx, f1(fx), linewidth=4)
plt.plot(fx, f2(fx), linewidth=4)
#Se agregan las leyendas de las curvas
plt.legend(["d=1","d=2"], loc="upper left")
plt.autoscale(tight=True)
plt.grid()
#Se muestra la gráfica en la pantalla
plt.show()

La gráfica que se obtienen es la siguiente:

figure_4

Si uno revisa el último código, se observa que se importan las librerías usuales, SciPy y Matplotlib, pero además la librería os (sistema operativo) y un módulo de nombre utils. Existe un módulo en Python con el nombre de utils, pero lo que realmente se está haciendo es usar un script escrito por Willi Richert y Luis Pedro Coelho, que son los autores de dicho libro, para extraer los datos de un archivo con el nombre de “web_traffic.tsv” desde algún directorio. Si solo se descargan los datos y se quiere correr el código en algún directorio  arbitrario donde no se encuentre ni los datos ni el archivo utils, Python no enviará un error.

Entonces para este último ejemplo es recomendable descargar el código en zip de GitHub del libro y revisar el código de los autores del libro.

Con estos ejemplos espero se de una idea de como se hace uso de Numpy, SciPy y Matplotlib. Sobre las otras librerías, en cada entrada comento lo necesario, pero cada uno tienen sus detalles al momento de usarlas. En algunos casos, como IPython, uno puede hacer uso de Python sin recurrir a su uso, pero no esta de más conocer el por qué es altamente usada y ya con la práctica terminaran convencidos de su utilidad.

Referencias:

1.-http://www.creativebloq.com/design-tools/data-visualization-712402

Libros:

1.-http://www.amazon.es/Building-Machine-Learning-Systems-Python/dp/1782161406

2.-http://www.amazon.es/Learning-Python-Mark-Lutz/dp/1449355730/ref=sr_1_1?s=foreign-books&ie=UTF8&qid=1431492984&sr=1-1&keywords=Python

3.-http://www.amazon.es/Python-Data-Analysis-Wrangling-IPython/dp/1449319793/ref=sr_1_6?s=foreign-books&ie=UTF8&qid=1431493009&sr=1-6&keywords=Python

4.-http://www.amazon.es/High-Performance-Python-Performant-Programming/dp/1449361595/ref=sr_1_11?s=foreign-books&ie=UTF8&qid=1431493009&sr=1-11&keywords=Python

5.-http://www.amazon.es/Natural-Language-Processing-Python-Steven/dp/0596516495/ref=sr_1_94?s=foreign-books&ie=UTF8&qid=1431493443&sr=1-94&keywords=Python

6.-http://www.amazon.es/Building-Machine-Learning-Systems-Python-ebook/dp/B00VAG2WU4/ref=sr_1_126?s=foreign-books&ie=UTF8&qid=1431493491&sr=1-126&keywords=Python