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

 

Comparación entre Máquina de Soporte Vectorial, Naive Bayes, Árboles de Decisión y Métodos Lineales

De lo que trata esta entrada.

En la entrada no explico a detalles técnicas sobre cada una de los algoritmos empleados, pero pueden consultarse en las referencias o en las categorías Machine Learning en R project y Machine Learning en Python.

Para el ejemplo solo uso código en R, la intención es usar tres muestras de datos simulados para mostrar las estimaciones que realiza la función svm de la librería e1071, solo en el caso de usarla para clasificar. Y por último uso unos datos  de correos clasificados como spam o no-spam provenientes del repositorio UCI Machine Learning Repository los cuales se encuentran cargados en la librería kernlab de  R project para aplicar las técnicas y comparar svm con otras técnicas.

Ejemplo.-Solo Suppor Vector Machine

En este ejemplo solo hago dos muestras de datos simulados, con distribuciones  Gaussianas y la intención es compartir los detalles de qué es en rasgos generales lo que hace el algoritmo de maquina de soporte vectorial al clasificar datos, esto es para dos dimensiones. Es decir, sobre un plano ya que visualmente se vuelve claro y con esa idea se puede pensar en el tipo de cosas que hace el algoritmo en más dimensiones.

Los primeros datos no se mezclan, por lo cual es visualmente claro como separarlos.

#Datos
#Librerías requeridas para el ejemplo

library(ggplot2)
library(e1071)

#Datos de mezcla 1

N<-500
x1<- rnorm(N) * 0.1
y1<- rnorm(N) * 0.1
X1<- t(as.matrix(rbind(x1, y1)))
x2<- rnorm(N) * 0.1 + 0.5
y2<- rnorm(N) * 0.1 + 0.5
X2<- t(as.matrix(rbind(x2, y2)))
X<- as.matrix(rbind(X1, X2))
X=data.frame(X)
X['Clus']=1
X[500:1000,3]=2
ggplot(X,aes(x=x1,y=y1))+geom_point(size=3,aes(colour=factor(Clus)))+ggtitle('Mezcla 1')+theme(plot.title=element_text(lineheight = 2,face='bold'))+
 xlab('Variable X1')+ylab('Variable Y1')

La gráfica de los puntos se ve así:Mezcla1_Gas

Para la siguiente mezcla considero tres muestras donde las clasifico en dos categorías.

#Mezcla 2
N<-500
x1<- rnorm(500,mean=0,sd=3)*0.1 
y1<- rnorm(500,mean=0,sd=3)*0.1 
X1<- t(as.matrix(rbind(x1, y1)))
x3<- rnorm(250,mean=0,sd=3)*0.1 
y3<- rnorm(250,mean=0,sd=3)*0.1+0.7 
X3<- t(as.matrix(rbind(x3, y3)))
x2<- rnorm(N,mean=1.5,sd=4)*0.1 + 0.5
y2<- rnorm(N,mean=0,sd=3)*0.1 + 0.5
X2<- t(as.matrix(rbind(x2, y2)))
Datos2<- as.matrix(rbind(X1, X3,X2))
Datos2=data.frame(Datos2)
Datos2['Clus']=1
Datos2[600:1250,3]=2
ggplot(Datos2,aes(x=x1,y=y1,colour=factor(Clus)))+geom_point(size=3)+ggtitle('Mezcla 2')+theme(plot.title=element_text(lineheight = 2,face='bold'))+
 xlab('Variable X1')+ylab('Variable Y1')

Mez2_GaussLa gráfica de los datos se ve como la anterior imagen.

Los datos muestran comportamientos distintos, los primeros es claro que se puede separar por una recta, pero la gráfica de la segunda muestra no se ve claramente que lo mejor sea separar los datos por una recta.

Como una buena práctica se debe de tomar una muestra de datos para entrenar el algoritmo (train set) y una cantidad de datos de prueba(test set). Para los siguientes ejemplos considero todos los datos, para ilustrar cual sería la gráfica de los datos predichos por el modelo.

La técnica de SVM para clasificar tienen de fondo la idea de encontrar el mejor “hiperplano” por medio del cual separar los datos. En este ejemplo el concepto de hiperplano es una recta, es decir; los datos están en dos dimensiones (x,y) y se busca la mejor recta que separe los datos, si uno piensa en datos con tres dimensiones (x,y,z) que pueden pensarse como alto, largo y ancho lo que se busca es encontrar el mejor plano que separa los datos, por ello el nombre de hiperplano.

Uso la librería e1071 para estimar el algoritmo SVM con el kernel lineal y radial para la primera muestra de datos, para la segunda uso lineal, polinomial, radial y sigmoid. La intención de este ultimo es compara cual de las cuatro muestra una similitud gráfica más parecida a los datos originales y comparo la tasa de predicciones correctas. Existen más librerías para hacer uso de SVM, dejo en la referencia las ligas.

Primera muestra de datos.

#SVM con Kernel Lineal
#Se procede a dejar una semilla para hacer una elección del mejor valor de parámetro cost

set.seed (1)
#Se hace cross-validation de k-fold, con k=10. Esto mediante la función tune()

tune.out=tune(svm,factor(Clus)~.,data=Datos1,kernel="linear",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)

#Se elige el mejor modelo con el mejor valor para el parámetro cost

bestmod=tune.out$best.model
summary(bestmod)
# Gráfica del mejor modelo con kernel lineal

plot(bestmod,Datos1)

#Predicción
Pred=predict(bestmod,Datos1)

table(predicción=Pred,Valores_reales=Datos1$Clus)

# Valores_reales
# predicción 1 2
#         1 499 1
#         2 0 500
#Con los valores obtenidos se tiene un 99.9% de eficiencia 

Algunas observaciones, el modelo pasa por un proceso de validación cruzada usando la función tune(), esto siempre es recomendable para elegir un buen modelo. El parámetro cost, significa el nivel de penalización que permite el modelo, en otras palabras el algoritmo busca la mejor recta que separe los datos y como tal el costo para encontrar esa recta requiere tolerar posibles datos que afectan a la estimación de dicha recta. Esto quizás es una explicación burda y mala, pero detrás de todo algoritmo de Machine Learning está un proceso de optimización el cual determina el valor de los parámetros que hacen que el algoritmo tenga el mejor valor posible.

La gráfica que se obtiene de este modelo es la siguiente:

SVM_Muestra1_ker-lineal

 

Esta gráfica muestra las dos clases que se buscaban definir.La eficiencia es muy alta, de 1000 datos clasifica correctamente el 99.9%.

Lo que se espera de otro kernel es que prevalezca la eficiencia de SVM y más aún que la gráfica de la clasificación sea casi igual a la obtenido cuando se usa un kernel lineal.

#Kernel radial
set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos1,kernel="radial",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_r=tune.out$best.model
summary(bestmod_r)
plot(bestmod_r,Datos1)
Pred=predict(bestmod_r,Datos1)
table(predicción=Pred,Valores_reales=Datos1$Clus)
# Valores_reales
# predicción 1 2
# 1 499 1
# 2 0 500
#Eficiencia de la clasificación 99.9%

Al obtener la gráfica del modelo por SVM se tiene:

SVM_Muestra1_ker-radial

Lo cual muestra una imagen muy similar a la obtenida con el kernel lineal. La tasa de eficiencia al clasificar es prácticamente la misma, 99.9%

Con la segunda muestra de datos pensar en separar por una recta no parece lo natural. Primero hago la estimación con cuatro tipo de kernels de la función svm y muestro la gráfica que regresa el modelo.

#Segunda muestra de datos
#Se usa primero el kernel lineal

set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="linear",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_l=tune.out$best.model
summary(bestmod_l)
plot(bestmod_l,Datos2)
Pred=predict(bestmod_l,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
# 1 504 99
# 2 95 552
#Eficiencia de la clasificación 84.4%

Se tiene la siguiente gráfica como resultado de implementar el algoritmo:

SVM_Mezcla2_ker-lineal

Si se compara los datos que predice el modelo con respecto a los originales se tienen lo siguiente:

Orig_vs_Ker-lineal

Se observa que la predicción muestra totalmente separadas las dos clases y cabe notar que se tiene el 84.4% eficiencia.

#Kernel polinomial

set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="polynomial",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_p=tune.out$best.model
summary(bestmod_p)
plot(bestmod_p,Datos2)
Pred=predict(bestmod_p,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
# 0 543 174
# 1 56 477
#Eficiencia de la clasificación 81.6%, con kernel polinomial

La gráfica que se obtiene del modelo es:

SVM_Mezcla2_ker-polin

Y otra vez haciendo una gráfica para comparar la predicción con los datos originales se tiene:

Orig_vs_Ker-pol

Se aprecia que no es tan definida la recta que separa las clases como en el ejemplo del kernel lineal, pero la eficiencia se reduce ya que es de 81.6%

#Kernel sigmoid

set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="sigmoid",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_s=tune.out$best.model
summary(bestmod_s)
#Informa del tipo de modelo, del que puede considerarse como el mejor modelo
plot(bestmod_s,Datos2)
Pred=predict(bestmod_s,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
#         1 509 102
#         2 90 549
#Eficiencia de la clasificación 84.6%, con kernel Sigmoid

La gráfica que se obtiene del modelo es la siguiente:

SVM_Mezcla2_ker-sigmoid

Haciendo la gráfica comparativa de predicción y datos originales se obtiene:Orig_vs_Ker-sigmoid

 

Se aprecia que es muy similar a la que se obtiene con el kernel lineal y hasta el nivel de eficiencia resulta muy aproximado, ya que es del 84.6%

#Kernel Radial
set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="radial",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_r=tune.out$best.model
summary(bestmod_r)
plot(bestmod_r,Datos2)
Pred=predict(bestmod_r,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
# 1 497 94
# 2 102 557
#Eficiencia de la clasificación 84.3%, con kernel lineal

La gráfica obtenida con este kernel es la siguiente:

SVM_Mezcla2_ker-radial

 

Y al comparar las predicciones con los datos originales se tiene:

Orig_vs_Ker-radial

Se aprecia que es similar a la obtenida por el kernel lineal y sigmoid, más aún la eficiencia resulta ser del 84.3%.

Entonces en resumen los datos al ser clasificados mediante SVM con distintos kernel resultó en estos datos resulta tener mejor eficiencia el kernel sigmoid, por un porcentaje mínimo sobre el lineal.

Para mostrar como se comporta SVM con una muestra de datos altamente no lineales o que resulta difícil separar por medio de una recta , genero una muestra más.

#Datos altamente no lineales

#Distribución uniforme 1250 valores
x1=runif(1250)-0.5
x2=runif (1250)-0.5
#Variable Indicadora
y=1*(x1^2-x2^2> 0)

#Gráfica donde se muestran los puntos pintados por etiqueta de y
#Se contruye un data.frame con los datos
X=data.frame(x1,x2,y)
ggplot(data=X,aes(x=x1,y=x2))+geom_point(aes(colour=factor(y),shape=factor(y)))+ggtitle('Datos Originales')+theme(plot.title=element_text(lineheight = 2,face='bold'))+
 xlab('Variable X1')+ylab('Variable X2')

La gráfica de los datos es la siguiente:

Datos_No-lineales

En esta muestra de datos se aprecia que separar las clases por una recta no parece lo más eficiente, para constatarlo estimo la clasificación con el kernel lineal.

#Kernel lineal
set.seed (1)
tune.out=tune(svm,factor(y)~.,data=X,kernel="linear",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod=tune.out$best.model
summary(bestmod)
#Informa del tipo de modelo, del que puede considerarse como el mejor modelo
plot(bestmod,X)
Pred=predict(bestmod,X)
table(predicción=Pred,Valores_reales=X$y)
# Valores_reales
# predicción 0 1
# 0 305 201
# 1 301 443
#Eficiencia de la clasificación 59.8%

La gráfica que se obtiene del SVM con kernel lineal es:

SVM_Mezcla3_ker-lineal

Haciendo la comparación entre datos originales y las predicciones se tiene lo siguiente:

Orig3_vs_Ker-lineal

Se aprecia que no es muy favorable usar el kernel lineal con estos datos, la eficiencia de la clasificación es del 59.4%

Para comparar con otro kernel uso el radial para hacer la clasificación.

#Kernel Radial

set.seed(1)
tune.out=tune(svm,factor(y)~.,data=X,kernel="radial",sacale=FALSE,type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod=tune.out$best.model
summary(bestmod)
#Informa del tipo de modelo, del que puede considerarse como el mejor modelo
Pred=predict(bestmod,X)
table(predicción=Pred,Valores_reales=X$y)
# Valores_reales
# predicción 0 1
# 0 603 7
# 1 3 637
#Eficiencia de la clasificación 99.2%, con kernel radial

La gráfica que se obtiene con este kernel es la siguiente:

SVM_Mezcla3_ker-radial

Esta imagen muestra que las curvas que se obtienen con el kernel radial son parecidas a las de los datos originales. Haciendo la comparación entre los datos y las estimaciones se obtiene lo siguiente:

Orig3_vs_Ker-Radial

Se aprecia que es muy similar y más aún la eficiencia es del 99.2%, lo cual comparado con el lineal es sumamente mejor usar el kernel radial.

Ejemplo con datos de correos clasificados

Los datos se encuentran cargados en la librería kernelabs, con la cual también se puede estimar la máquina de soporte vectorial. Las comparaciones las hago con varias técnicas:

  • Naive Bayes
  • LDA y Regresión Logística
  • Redes Neuronales
  • Árboles de decisión y random forest

Antes de eso aplico con los cuatro kernel disponibles en svm de e1071 para elegir el que mejor clasificación realiza en la muestra de prueba (test set).

Primero reviso los datos.

#Exploración breve de los datos
#Se cargan primero las librerías

library(kernlab)
library(e1071)
library(MASS)
library(nnet)
librery(randomForest)
data(spam)

#Datos
head(spam)
dim(spam)
#4601 58
str(spam)

#############################################
#Preparación de datos

index=sort(sample.int(4601,1150))
spam.train=spam[-index,]
spam.test=spam[index,]

Se puede hacer algo más al explorar la información, diseñar algunas gráficas o usar alguna técnica de estadística de exploración.

Lo que se hace en el código anterior es dividir los datos entre una conjunto de entrenamiento y un conjunto de prueba.

Aplico los cuatro kernel para elegir el que mejor clasificación realiza.

#SVM con diferentes kernel
########################################
#Kernel lineal

svmfit=svm(type~.,data=spam.train,kernel="linear",cost=0.1,type='C-classification' )
summary(svmfit)
Pred=predict(svmfit,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#         Valores_reales
#predicción nonspam spam
#   nonspam 646 40
#   spam     40 424
#Eficiencia del modelo 93.04%
#######################################
#Kernel radial
svmfit=svm(type~.,data=spam.train,kernel="radial",cost=0.1,type='C-classification' )
summary(svmfit)

set.seed(1)
tune.out=tune(svm,type~.,data=spam.train,kernel="radial",type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,10,100)))
tune.out
summary(tune.out)
bestmod_r=tune.out$best.model
summary(bestmod_r)
Pred=predict(bestmod_r,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#           Valores_reales
#predicción nonspam spam
#   nonspam     652 36
#   spam         34 428
#Eficiencia del modelo 93.91%
#################################################
#Kernel sigmoid
svmfit=svm(type~.,data=spam.train,kernel="sigmoid",cost=0.1,type='C-classification' )
summary(svmfit)
set.seed(1)
tune.out=tune(svm,type~.,data=spam.train,kernel="sigmoid",type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,10,100)))
tune.out
summary(tune.out)
bestmod_s=tune.out$best.model
summary(bestmod_s)
Pred=predict(bestmod_s,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#           Valores_reales
#predicción nonspam spam
#    nonspam    649 78
#    spam        37 386
#Eficiencia del modelo 90%
################################################
#Kernel polynomial
svmfit=svm(type~.,data=spam.train,kernel="polynomial",cost=0.1,type='C-classification' )
summary(svmfit)
set.seed(1)
tune.out=tune(svm,type~.,data=spam.train,kernel="polynomial",type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,10,100)))
tune.out
summary(tune.out)
bestmod_p=tune.out$best.model
summary(bestmod_p)
Pred=predict(bestmod_p,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#           Valores_reales
#predicción nonspam spam
#   nonspam     661 53
#   spam         25 411
#Eficiencia del modelo 93.2%

Se observa en las tasas de error que el mejor kernel de svm es el radial, así que considero este contra el cual comparo las otras técnicas.

Si bien no explico a detalle lo que hace el código en breve es estimar vía tune() la elección del mejor valor de cross-validation de varios valores.

Una técnica clásica para clasificar correos es Naive Bayes, la cual considera ciertas propiedades teóricas que permiten hacer un calculo rápido. Uso NaiveBayes de la librería e1071.

#Naive Bayes vs SVM elegido
NB=naiveBayes(type~.,data=spam.train)
Pred=predict(NB,spam.test)

table(predicción=Pred,Valores_reales=spam.test$type)
# Valores_reales
#predicción nonspam spam
# nonspam 362 18
# spam 324 446
#Se tiene una eficiencia del 70.2% de bien clasificados

Resulta que es mejor el modelo SVM ante el Naive Bayes por un porcentaje altamente mejor.

#Regresión Logística-LDA vs SVM elegido

Mod.lm=glm(type~.,data=spam.train,family='binomial')
summary(Mod.lm)
Mod.lmprob=predict(Mod.lm,newdata=spam.test,type="response")
Mod.pred=rep('nonspam',1150)
Mod.pred[Mod.lmprob>0.5]='spam'
table(Predicción=Mod.pred,Valores_Reales=spam.test$type)
#           Valores_Reales
#Predicción nonspam spam
#   nonspam 643 39
#   spam     43 425
#Predicción 92.86%
############################################
#Método LDA
lda.fit=lda(type~.,data=spam.train)
lda.pred=predict(lda.fit,spam.test)
lda.class=lda.pred$class
table(Predicción=lda.class,Valores_Reales=spam.test$type)
#lda.class nonspam spam
#   nonspam    650 91
#   spam        36 373
#Tasa de eficiencia 88.9%

Se tienen que el porcentaje de eficiencia en la clasificación de SVM 93.91% es mayor de casi 1% mejor que los métodos lineales, en este caso la regresión logística resulta ser la mejor entre las técnicas lineales.

Sin explicar detalles elijo un modelo de redes neuronales para clasificar, la librería es nnet y la función tienen el mismo nombre.

#Redes Neuronales
red=nnet(type~.,data=spam.train,size=2)
Pred=predict(red,spam.test,type="class")

table(Predicción=Pred,Valores_Reales=spam.test$type)
# Valores_Reales
#Predicción nonspam spam
# nonspam 650 29
# spam 36 435
#Tasa de eficiencia 94.34%

Resulta que el modelo de redes neuronales resulta ser mejor casi 0.5% mejor que el modelo de SVM. Lo cual faltaría hacer una revisión para definir el mejor modelo de redes neuronales y quizás la mejora es mejor o nula pero tendrías los valores óptimos de los parámetros.

#Dos técnicas de arboles
# La primera es por árboles y se hace cross-validation para elegir el mejor modelo
 
tree.spam=tree(type~.,spam.train)
set.seed (3)
cv =cv.tree(tree.spam,FUN=prune.misclass )

#Gráfica de Cross-validation para elegir el mejor valor de del parámetro best
par(mfrow =c(1,2))
plot(cv$size,cv$dev,type="b",col="2")
plot(cv$k,cv$dev,type="b",col="3")

prune.tree=prune.misclass(tree.spam,best=8)
tree.pred=predict(prune.tree,spam.test,type="class")
table(tree.pred,Valores_Reales=spam.test$type)

#         Valores_Reales
#tree.pred nonspam spam
#  nonspam 644 58
#  spam     42 406
#Tasa de clasificación 91.3%
#################################################
#Modelo de Random Forest

set.seed (1)
rf.spam =randomForest(type~.,data=spam.train,mtry=7, importance =TRUE)
rf.pred_spam=predict(rf.spam,newdata=spam.test,type="class")

table(rf.pred_spam,Valores_Reales=spam.test$type)
#            Valores_Reales
#rf.pred_spam nonspam spam
#     nonspam 663 29
#     spam     23 435
#Tasa de clasificación 95.4%

Revisando la tasa de clasificaciones correcta por medio de las técnicas de árboles de decisión , resulta ser la mejor técnica para clasificar los datos de los correos. Se obtienen que la técnica de Random Forest resulta clasificar correctamente el 95.4%, comparado con la red neuronal y SVM es mejor casi en un 2%.

La única intensión de este ejemplo es hacer una comparación entre diversas técnicas de clasificación, de cierto modo hacer un ejemplo sencillo y que se puede replicar. Ha esto le faltaría agregar algunos gráficos que muestran los resultados o hacer algunas variaciones de los parámetros requeridos en cada algoritmos para determinar con mejor eficiencia el mejor modelo.

Espero ilustre de manera breve la variedad de técnicas y permita visualizar que al final se requiere probar con muchos modelos para elegir uno entre todos los realizados.

Algoritmos de Machine Learning en R project

“The best thing about R is  that it was developed by statisticians. The worst thing about R is that….it was developed by statisticians.”  Bo Cowgill, Google, Inc.

Todas las entradas de esta categoría muestran algoritmos y técnicas clasificadas como algoritmos de Machine Learning o Aprendizaje Estadístico. En cada entrada comparto el código  y algunos repositorios donde pueden bajar datos y notas sobre los temas.

Las ventajas de usar R project para dar estos ejemplos son las siguientes:

1.-R project es software libre u open source, por lo cual cualquier persona pude descargarlo.

2.-El código que comparto solo debe de ser copiado y pegado en un Script o en la consola  de R project para reproducir los resultados.

Otra ventaja de aprender estas técnicas con ejemplos en R, es que muchos de los algoritmos son acompañados de alguna técnica gráfica para visualizar el resultado. R project por default cuenta con una gran cantidad de herramientas gráficas o se puede hacer uso de librerías como ggplot2, lattice, ggvis para hacer gráficas aún más elaboradas y visualmente más agradables.

El entendimiento de la técnica después de realizar un ejemplo desde mi perspectiva genera interés por saber como aplicarla a otros datos y nos invita a conocer más sobre el algoritmo desde un punto de vista no solo práctico, sino también teórico.

Evité en la medida de lo posible hacer uso de ecuaciones y en cada ocasión que se menciona algún concepto de matemáticas trato de explicarlo de un modo sencillo. Si este concepto requiere a mi parecer una mejor lectura, agregue algunas ligas para consultar en la red.

La única meta que me propuse fue explicar las cosas de modo tal que cualquier persona con un poco de curiosidad pueda entender algo de cada entrada.

Teniendo en mente lo anterior, espero que  el objetivo de trasmitir el conocimiento desde un ejemplo concreto haya sido logrado. En caso de que no sea así, sería grato saber qué parte no es clara y como mejorar este material.

Resumen de las entradas.

El orden siguiente es el que considero apropiado para revisar cada entrada.

1.-¿Qué se necesita conocer de R project para aprender algunas técnicas de Machine Learning?

Es una breve explicación sobre R project y sobre algunas de las funciones con las que cuenta por default. También menciono las librería que se requieren para replicar los ejemplos, aunque en cada entrada menciono algo breve sobre la librerías requeridas.

2.-Análisis Confirmatorio vs Exploratorio

Explico la diferencia entre los dos tipos de análisis y cómo están relacionados. Clasifico el tipo de medidas que uno puede analizar en los datos y pongo algunos ejemplos del tipo de gráficas que pueden ayudarnos a visualizar las medidas al momento de explorar la información.

3.-Regresión Lineal, un ejemplo sencillo.

Esta técnica es muy sencilla de implementar y los datos que nos regresa el comando en R debe ser entendidos para poder evaluar la calidad del resultado obtenido. Omito detalles formales del tema y comparto algunos ejemplos de regresiones simples y múltiples. Es un tema con muchos detalles importantes, por ejemplo la revisión de las posibles fallas al construir un modelo lineal o multilineal. No toco mucho  esos temas, pero en las referencias pueden encontrar observaciones y técnicas para conocer más al respecto.

4.-Regresión no lineal, Cross-Validation y Regularization. 

Los tres temas pueden haber requerido una entrada para cada uno, pero funciona bien mostrar como se relacionan. El objetivo es mostrar la implementación de la regresión no lineal y validar que el modelo elegido es el adecuado, para evitar la sobre estimación de los datos. En esto último las cross-validation y la regularización ayudan a tener un método para hacer una elección del modelo. Se dan algunos ejemplos del uso de los tres conceptos y técnicas.

5.-Clasificación binaria….Naive Bayes. 

Es uno de los algoritmos más usados por sencillo y por que puede ser modificado para tener una clasificación rápida en varias categorías. Doy el ejemplo clásico sobre detección de Spam y trato de complementarlo con el uso de una librería en la cual se cuenta con una función para estimar el modelo. Agregué una pequeña explicación sobre la probabilidad condicional y sobre el uso de la librería tm para analizar textos en R project. Espero ser claro y que permita comprenderse el ejemplo.

6.-Análisis de Cluster, un ejemplo sencillo. 

Se explica de modo breve como funciona la detección de Clusters. Como el tema requiere hacer uso de la noción de distancia o métrica, trato de explicar el concepto. Omití mucho material, sobre todo en el uso de distintas métricas para distintos tipos de variables.

7.-Un ejemplo sencillo sobre Análisis de Componentes Principales (PCA, Principal Component Analysis).

La técnica es muy usada y conocida para hacer reducción de dimensiones, explico qué es reducir dimensiones y como interpretar los resultado de la técnica. Como ejemplo principal hago una reducción de 25 indicadores para formar uno que resuma la información del resto. Replico un ejemplo de análisis de la digitalización de los números escritos a mano, este es un ejemplo clásico y la base de datos fue objeto de estudio para comparar técnicas de clasificación.

8.-Sobre sistemas de recomendación, ejemplo sencillo. 

Hago un breve ejemplos de un sistema de recomendación usando la técnica de k-vecinos cercanos (K-nn), esta técnica es no-paramétrica. Para ejemplificar como funciona pongo un ejemplo inicial y trato de explicar qué realiza el algoritmo y aplicarlo a un caso concreto. Los sistemas de recomendación prácticamente dan la libertad de elegir la técnica de clasificación que uno considere adecuada para el tipo de datos o servicio, por lo cual existen libertad de replicar este ejemplo haciendo uso de otro técnicas de clasificación.

9.-Optimización.

Esta entrada creo que es extensa y puede ser un poco amplia,ya que traté de explicar desde aspectos básicos de mínimo y máximos hasta algoritmo estocásticos y genéticos. Puse un ejemplo de como usar el código y en los casos más sofisticados, la parte estocástica, comparto un ejemplo que permitiera comparar los resultado obtenidos con diversos tipos de algoritmos, no puse los tiempos que tardan en dar la solución pero en caso de replicarlos será notorio el costo del procesamiento.

10.-Máquina de soporte Vectorial (SVM-Sopport Vector Machine).

La técnica sumamente usada y funciona bien para datos no-lineales, se dan varios ejemplos de como usar dos librerías que permiten calcular SVM con distintos kernels. En uno de los ejemplos de hace la comparación de varios modelos para determinar por medio de MSE cual es el mejor. Son ejemplos muy sencillos y trato de agregar explicaciones a cada ejemplo.

Lo que no está aún en estas entradas, pero que estará en inicio de Marzo 2016.

Dejé fuera por el momento varios temas fuera de las entradas, ejemplo Redes neuronales , Arboles de decisión, Random Forests, AdaBoost. Espero posteriormente agregar entradas sobre esos algoritmos.

En general todos las entradas tratan de simplificar y resumir lo que se puede encontrar en el texto  “Machine Learning for Hackers“, libro escrito por Jhon M. White y Drew Conway. Si bien, el libro en si es sencillo y requiere un bajo nivel de matemáticas para estudiarse, traté de hacer más sencillas las explicaciones y de complementarlo con otros ejemplos.

La recomendación es complementar las explicaciones de las entradas con la lectura del texto.  Jhon M White compartió su código  y datos usados en los ejemplos en el repositorio Github y pueden descargar el código presionando estas letras en rojo.

Algunos de los códigos deben de revisarse ya que las librerías que usaron en el año de la edición del texto ya fueron actualizadas, por lo cual es generan algunos errores.

Dos temas que vienen en el texto y decidí omitir son generación de un rango y análisis de redes sociales o de grafos. En otro momento comparto ejemplos respecto al tema.

El segundo es sumamente interesante; análisis de redes sociales, pero el código con el que cuenta el texto ya no es funcional. Así que se requiere rehacer buena parte del código y diseñar nuevos ejemplos, para eso se necesita cierto tiempo. Más aún debido al boom de las redes sociales existen otro tipo de herramientas sobre las cuales se puede escribir.

Para tener una idea de lo nuevo e interesante se pueden consultar las siguientes ligas para echar un vistazo:

Deseo que el material que se encuentra en esta categoría les sea útil y sobre todo les deje con ganas de aprender más al respecto y no solo desde el lado aplicado, sino también desde el lado teórico.

En la siguiente cita aplique “instrumento= algoritmo”.

“Tienes que aprender a tocar tu instrumento. Después debes practicar, practicar y practicar. Y después, cuando finalmente estás en el escenario, olvídalo todo y ulula.” Charlie Parker

 Nota: Todas las críticas y comentarios respecto a las entradas son bienvenidos.