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

 

Optimización

Sobre Optimización

Para entender en general las ideas de esta entrada basta con tener algunos conocimientos de Calculo Diferencial, de una y varias variables. Si no se cuenta con ideas o conocimientos de ello, no implica el no entender los algoritmos de optimización que se comentan o por lo menos visualizar como funcionan.

Parto de ideas sencillas hasta dar unos ejemplos de optimización estocástica. Espero poder trasmitir la idea de los métodos y que puedan servirles de ayuda para investigar más respecto al tema.

Omito hablar sobre el proceso ” gradiente descendiente estocástico” pero dejo referencias donde pueden consultar en qué consiste este método y por que es que se hace amplio uso en técnicas de Machine Learning.

 La idea de optimizar  seguro es familiar cuando vemos como nuestras madres resuelven la siguiente pregunta: ¿cómo hacer para comer carne con solo cierta cantidad de dinero?

Aunque parece broma la pregunta anterior es una optimización, tenemos que encontrar las cantidades de bienes que podemos comprar considerando que solo tenemos una cantidad de dinero y debemos de incluir en la lista de bienes la carne. Obviamente esto no dice que nuestras madres hagan uso de las técnicas que comento esta entrada, pero espero que de una idea burda de lo que es optimizar.

Otro ejemplo es cómo maximizar las ganancias de una empresa o negocio, considerando R cantidad de ingresos y C cantidad de costos. Claramente este problema implica conocer qué genera nuestros ingresos y costos, para plantear un modelo idealizado sobre la empresa.

En general la idea de optimizar es encontrar los máximos o mínimos de una función. Para tener una idea de qué significa hago dos gráficas como ejemplo en R project.

#Creamos una función
#Ejemplo 1-Mínimo de una función

Func_1<-function(x){
 x[1]^2 
      }

resultado<-optimize(Func_1,c(-5,5),maximum=FALSE,tol=1e-8) 

#Vemos los resultados
print(resultado$minimum)
print(resultado$objective)

#Generamos el rango de valores para x y y=f(x)
x<-seq(-5,5,length.out=100)
y<-Func_1(expand.grid(x))

#Graficamos identificando el mínimo

plot(x,y,xlab="x",ylab="f(x)",col="4")
points(resultado$minimum,resultado$objective, col="red",pch=19)
rect(resultado$minimum-0.3, resultado$objective-0.7, resultado$minimum+0.3,resultado$objective+0.7)

Ejemplo_opt_mínimo

La gráfica muestra en color rojo el mínimo de la función, lo obtenemos por medio de la función optimize. La cual necesita una función objetivo, un rango de valores para optimizar y un nivel de tolerancia en el cálculo. Esto último es el número que utilizará para realizar su aproximación. Recomiendo revisar la información de la función por medio del comando help.

Cuando pedimos el resultado por medio de print(resultado$minimun)en el código, se observa que el valor es 0, lo cual es evidente al ver la gráfica. Hago otro ejemplo con una función diferente.

#Ejemplo 2
#Creamos una función
Func_2<-function(x){
 x^2+6*sin(x*pi*2)
 }
 
#Optimizamos
resultado2_1<-optimize(Func_2,c(-5,5),lower="-1",upper="3",maximum=FALSE,tol=.Machine$double.eps^0.25) 
resultado2<-optimize(Func_2,c(-5,5),maximum=FALSE,tol=1e-8)
#Vemos los resultados
print(resultado2$minimum)
print(resultado2$objetive)

#Rango de calores para nuestras gráficas 
x2<-seq(-5,5,length.out=100)
y2<-Func_2(expand.grid(x2)) 

#Graficamos
plot(x2,as.matrix(y2)[,1],xlab="x",ylab="f(x)",col="4")
points(resultado2$minimum,resultado2$objective, col="red",pch=19)
rect(resultado2$minimum-0.3, resultado2$objective-0.7, resultado2$minimum+0.3,resultado2$objective+0.7)

points(resultado2_1$minimum,resultado2_1$objective,col="green",pch=19)
rect(resultado2_1$minimum-0.3,resultado2_1$objective-0.7,resultado2_1$minimum+0.3,resultado2_1$objective+0.7)

Ejemplo_opt_mínimo_2

La diferencia de esta última gráfica es que se tienen dos puntos, rojo y verde, que a primera vista parece que el verde es el “mínimo” de toda la gráfica. Pero lo que realmente hice fue acotar el cálculo del optimización, por lo cual si observamos en la función  optimize hice distintas cotas (lower, upper) para que se ejecutara la función. Entonces lo que se ve son mínimos locales y dependen de la restricción que se haga a los cálculos.

¿Qué es lo que hace la función optimize? Da el punto máximo o mínimo según las restricción para el cálculo o en otro caso regresa el máximo o mínimo global.

La diferencia en local y global, es que el primero se restringe a una pequeña región, ejemplo a las x que están entre -2 y -1, el punto rojo de la gráfica es un mínimo local. Por otro lado, global se refiere a que para todo el rango de las x’s que se revisan no se puede encontrar otro valor de f(x) menor. En el ejemplo el punto verde de la segunda gráfica es un mínimo global y el punto rojo de la primera también.

Los dos ejemplos muestran gráficamente cual es la idea de obtener un mínimo que depende de una sola variable x. Pero en la mayoría de problemas de optimización no son funciones de una sola variable, sino de varias variables. Ejemplo, cuando  se busca maximizar ganancias o minimizar costos, es difícil que dependa de un solo producto que se vende o el costo de un solo insumo.

Entonces, la idea de optimizar se puede generalizar para varias variables. Pero en este caso se tienen varios métodos para hacer el cálculo y cada uno tienen ventajas y desventajas, depende del tipo de problema o función que  se desea optimizar.

¿Qué tiene que ver todo esto con Machine Learning?

 En general todos los algoritmos de Machine Learning tienen tres piezas, datos de entrada, datos de salida y una función entre los datos. Pero por medio de un conjunto de datos de entrenamiento y un algoritmo de aprendizaje se trata de obtener la mejor función entre los datos.

Entonces, la optimización está ligada a los algoritmos, ya que siempre en los algoritmos  se tiene una función objetivo y un conjunto de parámetros que se busca estimar. Por lo cual es fundamental la optimización al aplicar los algoritmos de Machine Learning.

Por ejemplo, en la regresión lineal. Sobre la cual en otra entrada se mencionó que la función lm permite estimarla. Se toman datos de altura y peso de personas, los datos se pueden descargar desde de Github.

Al usar la función lm  no se muestra lo que hace o como estima los valores de la regresión  a y b, donde la ecuación es:

Peso=a+b*Altura

Lo que se puede hacer es estimar los valores de los parámetros por medio de una optimización y comparar con los valores que se estiman en R project por medio de la función lm.

#Ejemplo de estimación de Parámetros
#Definimos la función de la regresión

Altura.a.peso <- function(altura, a, b)
{
 return(a + b * altura)
}

#Función de error cuadrado
error.cuadrado <- function(altura.peso, a, b)
{
 predic <- with(altura.peso, Altura.a.peso(Altura, a, b))
 errors <- with(altura.peso, Peso - predic)
 return(sum(errors ^ 2))
}

#Optimización

optim(c(0, 0),function (x){error.cuadrado(altura.peso, x[1], x[2])})

El código anterior hace dos cosas. La primera es definir la función de la cual se busca estimar sus parámetros a y b, pero esa función solo define rectas y para elegir la que mejor describe la relación entre  los datos se necesita otra función, error.cuadrado.

La función error.cuadrado como se observa depende de la función de Altura.a.peso. Así que se tiene solo una función compuesta con otra y  se busca en qué valores se minimiza esa función.

La función error.cuadrado permite encontrar la mejor recta entre los valores de Peso contra Altura. Además indica que se encontró la recta que tienen el menor error cuadrado medio, que es medido entre el Peso menos el valor que estimamos mediante la función Altura.a.peso.

Comparo los valores que se encuentran con la estimación por optimización y por la función lm.

#Comparación entre las estimaciones

 >lm(Peso~Altura,data=altura.peso)

Call:
lm(formula = Peso ~ Altura, data = altura.peso)

Coefficients:
(Intercept) Altura 
 -350.737 7.717 


> optim(c(0, 0),function (x){error.cuadrado(altura.peso, x[1], x[2])})
$par
[1] -350.786736 7.718158

$value
[1] 1492936

$counts
function gradient 
 111 NA 

$convergence
[1] 0

$message
NULL

Se nota que los valores de la regresión por medio de lm y los valores estimados por la función optim son prácticamente los mismos. Entonces se obtiene mediante optimización de una función objetivo error.cuadrado  el valor de sus parámetros a, b.

Se observa que la optimización fue hecha para una función de dos variables, el método que se usa para optimizar es Nelder-Mead, más adelante hablo de él con algo de detalle.

En conclusión, lo que parecía que operaba como una caja-negra cuando se usa la función lm, ahora se puede entender o tener una visión de cómo posiblemente es el proceso de estimación de los parámetros.

Aprovechando el ejemplo de la regresión lineal se puede modificar la función de error para estimar lo que se llama Rigde Regression. Es otro tipo de regresión en la cual para estimar los parámetros se afecta su error cuadrado o se penalizan los errores. El objetivo de esta técnica es evitar la sobre estimación de los datos y disminuir la varianza del ajuste lineal.

La Ridge Regression implica agregar un nuevo parámetro, lambda; el cual para hacer una buena elección se requiere hacer uso de la Regularización y Cross-Validation[1]. Se haciendo la modificación a la estimación para revisar el comportamiento al optimizar.

Haciendo dos estimaciones de Ridge Regression y se compara la estimación con respecto a las que se realizo anteriormente.

#Calculo de Ridge Regression

ridge.error <- function(altura.peso, a, b, lambda)
{
 predic <- with(altura.peso, Altura.a.peso(Altura, a, b))
 errores <- with(altura.peso, Peso - predic)
 return(sum(errores ^ 2) + lambda * (a ^ 2 + b ^ 2))
}

lambda=1
optim(c(0, 0),function (x){ridge.error(altura.peso, x[1], x[2], lambda)})

$par
[1] -340.434108 7.562524
lambda=0.75
optim(c(0, 0),function (x){ridge.error(altura.peso, x[1], x[2], lambda)})

$par
[1] -343.13351 7.60324

El ejemplo muestra el uso de la optimización al estimar los parámetros de un modelo. En la Rigde Regression se definió el valor de lambda, lo cual hace que se tenga ese parámetro “libre”  en el modelo y por ello es necesario para un buen uso de esta técnicas elegir la mejor lambda por medio de regularización y cross-validation.

Entonces al usar el paquete ridge de R project para estimar la ridge regression, ya no parece del todo extraño como hace el cálculo el software. Podría usar la función linearRidge del paquete ridge y comparar por medio de optimización, regularización y cross-validacón la elección del lambda.

Algo que no hice pero que puede resultar ser un buen ejercicio es revisar el artículo de investigación de Erika Cule y estudiar el método de elección de lambda contra otro método para elegirla.

Ejemplo_opt_regresion

La gráfica muestra como cambia ligeramente la pendiente de las rectas. Haciendo solamente una modificación en la estimación de los errores, regresión lineal y ridge regression.

Espero que este ejemplo haga notar que los algoritmos de Machine Learning  que parecen cajas negras, en el fondo están combinando cosas, y una de esas es un proceso de optimización para estimar los parámetros.

Este tema un trabajo serio al momento de procesar volúmenes masivos de datos, quienes tengan interés pueden revisar las investigaciones del grupo de optimización de la universidad Carnegie-Mellon.

Cabe mencionar que en las investigaciones sobre optimización su relevancia puede implicar una disminución en el costo computacional o en la mejora de algoritmos.

Ejemplos de técnicas de Optimización

Lo siguiente que comparto en esta entrada son técnicas de optimización convexa y dos técnicas de optimización estocástica. Para algunas solo pondré ejemplos de cómo se estima y como interpretar el resultado siguiendo unos ejemplos sencillos publicados por Jason Brownlee en su página. En otros pondré ejemplos de como usarlas de manera más detallada.

La explicación de cada uno de los métodos no es exhaustiva, tanto por que está fuera de mi alcance como de la intensión de esta entrada.

Dejo referencias en los comentarios y al final de la entrada algunas ligas a notas, reportes de investigación y libros para consultar más sobre el tema.

Método 1.-Búsqueda de la sección dorada

No repito el código, pero la primera gráfica que se hizo en esta entrada hacía uso de la función optimize y muestra un punto rojo que es el mínimo de la función. El método por medio del cual estima ese punto es conocido como  búsqueda de la sección dorada.

Su uso se limita a función que dependen de una sola variable y se busca el punto óptimo entre dos puntos que se definen del rango de las x’s. El nombre provienen de la relación entre los tres primeros puntos que toma el método, ya que esto tienen una proporción con el número conocido como razón dorada o número áureo.

Entre los limitantes de dicha técnica es que solo funciona para funciones de una sola variable, de preferencia que sea una función convexa y que se pueda identificar que se tienen un único óptimo entre la región que define para realizar el cálculo.

Función Rosenbrock

Para tres de los siguientes ejemplos se hace uso de la función Rosenbrock. Es una función no convexa que se usa para probar los métodos de optimización o rendimiento de los algoritmos. En general se puede usar de dos variables o se puede generalizar para n-variables. Forma parte de un conjunto de funciones que son usadas para hacer pruebas con técnicas de optimización.

La ecuación en general de la función para dos variables es:

f(x,y)=(a-x)^2+b(y-x^2)^2

La visualización de la función para un par de valores a=1 y b=100.

#Función Rosenbrock
#La función es f(x,y)=(1-x)^2+100*(y-x^2)^2

nx<-21
ny<-21

x<-seq(-5,5,length=nx)
y<-seq(-5,5,length=ny)

#Se genera la salida de los valores en z

z<-outer(x,y,function(x,y)(1-x)^2+100*(y-x^2)^2)

#Escala de colores
hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1,-1])
hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt))

#Se construye la gráfica en perspectiva con cierto ángulo para visualizarla
persp(x, y, z, col =cm.colors(10)[floor(9* hgt + 1)],theta = 35,phi=30)

Función_Rosenbrock

Método 2.-Nelder-Mead

El método es adecuado para funciones no lineal de varias variables y se busca el óptimo sin restricciones en el dominio. La herramienta que usa para buscar el óptimo en su algoritmo es una estructura geométrica conocida como simplex. El método para una función de dos variables parte de la elección de 3 puntos iniciales de los cuales bajo un proceso de varias etapas se va modificando el simplex hasta aproximarse al punto óptimo de la función. El proceso de transformación del simplex es la parte principal del algoritmo.

Es adecuado para funciones de varias variables y al no requerir la derivada para su algoritmo puede ser usada para funciones no diferenciables, no continuas o con ruido. Un resumen de sus ventajas y desventajas se pueden consultar en, Nelder-Mead_algorithm.

#Creamos una función

#Definición de una función Rosenbrock, el optimo de ubica en (1,1)

rosenbrock <- function(v){ 
 #Valores de a=1, b=100
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}

# Identificación del mínimo local por el método Nelder-Mead 
resultado3 <- optim(c(runif(1,-3,3), runif(1,-3,3)),rosenbrock,NULL,method="Nelder-Mead",control=c( maxit=100,reltol=1e-8,alpha=1.0, beta=0.5,gamma=2.0))

#Revisión de resultados
print(resultado3$par) 
print(resultado3$value) 
print(resultado3$counts)

#Gráfica de cortorno
x3<- seq(-3, 3, length.out=100)
y3<- seq(-3, 3, length.out=100)
z3<- rosenbrock(expand.grid(x3, y3))
contour(x3, y3, matrix(log10(z3), length(x3)), xlab="Valores de X",ylab="Valores de Y")
#Punto óptimo
points(resultado3$par[1], resultado3$par[2], col="red", pch=19)
#Úbicación del punto óptimo
rect(resultado3$par[1]-0.2, resultado3$par[2]-0.2, resultado3$par[1]+0.2, resultado3$par[2]+0.2, lwd=2)

Ejemplo_opt_Nelder-Mead2

Método 3.-Gradiente Descendente o Algoritmo de Descenso Rápido.

El métodos hace uso del gradiente de una función y cuando busca el mínimo hace uso del negativo del gradiente para aproximarse al óptimo.

Funciona bien si la función es convexa, ya que el método estimará el mínimo global. El método es muy sensible a la elección del punto inicial para hacer la estimación y como depende de un parámetro para ir haciendo los avances, es susceptible a malas elecciones de él.

Como referencias dejo las siguientes ligas: Aplicación del método a procesamiento de imágenes y curso de Andrew Ng en línea.

Usamos el método para calcular el mínimo de una función que tienen su mínimo global en el punto (0,0), el cual se visualiza en la imagen.

Ejemplo_función_min_global

#Creamos una función
#Se define una función, el óptimo se ubica en (0,0)

Func_3<- function(x) {
x[1]^2 + x[2]^2
}

#Definición de la derivadas
derivada <- function(x) {
c(2*x[1], 2*x[2])
}

#Definición del método del gradiente
gradient_descent <- function(func, derv, start, step=0.05, tol=1e-8) {
 pt1 <- start
 grdnt <- derv(pt1)
 pt2 <- c(pt1[1] - step*grdnt[1], pt1[2] - step*grdnt[2])
 
 while (abs(func(pt1)-func(pt2)) > tol) {
 pt1 <- pt2
 grdnt <- derv(pt1)
 pt2 <- c(pt1[1] - step*grdnt[1], pt1[2] - step*grdnt[2])
 print(func(pt2)) #regresa el valor en el punto
 }
 pt2 # Regresa el valor óptimizado
 }

#Método del Gradiente Descendente
resultados4 <- gradient_descent(Func_3, derivada,c(runif(1,-3,3), runif(1,-3,3)), 0.05,1e-8)

#Resumen de resultados
print(resultados4) 
print(Func_3(resultados4)) 

#Gráfica de contorno
x4 <- seq(-3, 3, length.out=100)
y4 <- seq(-3, 3, length.out=100)
z4 <- Func_3(expand.grid(x4, y4))
contour(x4, y4, matrix(z4, length(x4)), xlab="x",ylab="y")
#Dibujamos el punto óptimo
points(resultados4[1], resultados4[2], col="red", pch=19)
#Dibujamos el cuadrado del óptimo
rect(resultados4[1]-0.2, resultados4[2]-0.2, resultados4[1]+0.2, resultados4[2]+0.2, lwd=2)

Ejemplo_gradiente_descendente

Método 4.-Gradiente Conjugado

Es adecuado para funciones no lineales de varias variables. Le idea del método es construir una base ortogonal de vectores y hacer uso de ella para hacer la búsqueda del óptimo de un modo más eficiente. Aspectos teóricos de manera breve puede consultar en las notas en español del Prof. Maximenko y una buena referencia para ver más detalles recomiendo las notas del Prof.Shewchuk.

Igual que otros métodos su eficiencia se ve mermada ante una mala elección del punto inicial en  funciones no  convexas, pero una de sus ventajas es que no hace uso de la matriz Hessiana en su algoritmo lo cual permite hacer uso del método para funciones de muchas variables de manera eficiente.

Para hacer el ejemplo hacemos uso de la función Rosenbrock.

#Definimos la función Rosenbrock

rosenbrock <- function(v) { 
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}

#Definición de la derivada de la función
derivada2 <- function(v) {
c(-400 * v[1] * (v[2] - v[1]*v[1]) - 2 * (1 - v[1]), 
200 * (v[2] - v[1]*v[1]))
}


Ejemplo_opt_gradente_conjugado2

Método 5.-BFGS

Este método hace uso tanto del gradiente como de una aproximación a la inversa de la matriz Hessiana de la función, esto para hacer una aproximación al cálculo de la segunda derivada. Por ser un método de aproximación a la segunda derivada se dice que es un método Quasi-Newtoniano. El problema con el método es el costo computacional para funciones de muchas variables, ya que se almacena una matriz cuadrada de datos tan grande como el cuadrado de la cantidad de variables.

Es adecuado para funciones no lineales de varias variables y la búsqueda del óptimo sin restricciones. Si bien el método hace que la convergencia al óptimo sea rápida por ser de aproximación a la segunda derivada, el costo de procesamiento es bastante alto.

Se pueden consultas las notas de Galen Andrew o de Peter Blomger, para ver detalles del método y el algoritmo.

#Método BFGS

#Definición la función Rosenbrock, se óptimiza en (1,1)
rosenbrock <- function(v) { 
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}
 
#Definición del Gradiente
derivada4 <- function(v) {
c(-400 * v[1] * (v[2] - v[1]*v[1]) - 2 * (1 - v[1]), 200 * (v[2] - v[1]*v[1]))
}
 
#Calculo del método BFGS 
resultados6 <- optim(c(runif(1,-3,3), runif(1,-3,3)),rosenbrock, derivada4,method="BFGS",control=c(maxit=100,reltol=1e-8))
 
#Resumen de resultados
print(resultados6$par)
print(resultados6$value) 
print(resultados6$counts)
 
#Gráfica de contorno
x5 <- seq(-3, 3, length.out=100)
y5 <- seq(-3, 3, length.out=100)
z5 <- rosenbrock(expand.grid(x5, y5))
contour(x5, y5, matrix(log10(z5), length(x5)), xlab="x", ylab="y")

points(resultados6$par[1], resultados6$par[2], col="red", pch=19)
# draw a square around the optima to highlight it
rect(resultados6$par[1]-0.2, resultados6$par[2]-0.2, resultados6$par[1]+0.2, resultados6$par[2]+0.2, lwd=2)

Ejemplo_BFGS

Los métodos anteriores son considerados de optimización convexa, para ver detalles teóricos recomiendo consultar un texto especializado en el tema. Una buena referencia es el libro de Stephen P. Boyd y Lieven Vandenberghe el cual se puede descargar libremente desde su página. El segundo autor tienen dos cursos en línea con formato de presentaciones; optimización convexa y métodos de optimización para sistemas de larga escala.

Optimización Estocástica (o ¿computación evolutiva?)

Los siguientes ejemplos que comparto son abordados con técnicas de optimización estocástica, las cuales también son estudiadas en computación evolutiva. Por supuesto que no doy explicaciones exhaustivas, pero dejo referencias al final de la entrada.

Los métodos estocásticos que generalmente se usan son: algoritmos genéticos (genetic algorithm)  y recorridos simulados (Simulated anneling).

En el ejemplo considero el problema siguiente, si se  tiene un grupo de viajeros y todos parten de lugares distintos pero se dirigen al mismo sitio, si todos están de acuerdo en reunirse en un punto y se esperaran para compartir transporte terrestre, esto implica que entre los posibles vuelos tienen que elegir aquellos que mejor les convenga en costo, tiempo de espera, horarios de llegada y salida. Así que se tienen varios factores a considerar en la optimización.

Como cualquier problema de optimización se requiere tener la función objetivo, en este ejemplo se le llamará función de costo.

Los datos pueden ser descargados desde la página de Toby Segaran, cuenta con con una versión del código en python. Hice una adaptación en R y considerando las restricciones del problema y de los datos.

Los pasos para abordar el problema son los siguientes:

  1. Se representarán las soluciones como cadenas de números.
  2. Se presentaran las soluciones como un tabla.
  3. Se definirán una función objetivo o función de costo.
  4. Se estimaran soluciones por medio del método búsqueda aleatoria, recorrido simulados y algoritmos genéticos.

Paso 1

Se representa las solución como una cadena de 12 dígitos que pueden estar entre 1 y 9. Ejemplo, solución=[1,3,4,5,8,10,2,3,7,9,3,5], representa que la persona uno toma el primero vuelo y regresa en su tercer vuelo posible. La persona dos toma el cuarto vuelo y regresa en su quinto vuelo posible, así sucesivamente para las 6 personas.

#Carga de los datos desde el directorio
viajes.path<-file.path("data","schedule.txt")

viajes<-read.table(viajes.path,sep=",",stringsAsFactor=FALSE)

colnames(viajes)<-c('Origen','Destino','Hora de Salida','Hora de arrivo','Precio')

#Generamos lista de viajeros

gente=list(c('Seymour','Franny','Zooey','Walt','Buddy','Les'),c('BOS','DAL','CAK','MIA','ORD','OMA'))

#Destino al que viajan
Destino=c('LGA')

#Definimos una función para medir el tiempo, será necesaria en otras funciones
#Librería para manipulas fechas
library(lubridate)

#Función que nos regresa los minutos
minutos<-function(horario){
 min<-hour(hm(horario))*60+minute(hm(horario))
 return (min)}

Paso 2

Teniendo los datos cargados, se procede a representar por una pequeña tabla la solución al problema. Para ello se define la función impri.horario.

#Impresión de soluciones

impri.horarios<-function(r){

#r es la cadena con datos de los viajeros y lugars
 for( i in 1:(length(r)/2)){
 nombre=gente[[1]][i]
 origen=gente[[2]][i]
 Destino=c('LGA')
 
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 
 a<-c(nombre,origen,paste(salida[,3],"-",salida[,4],sep=""," ",salida[,5]),paste(regreso[,3],"-",regreso[,4],sep=""," ",regreso[,5]))
 print(a,sep="\t")
 
 } 
 
} 

Paso 3

Se define una función objetivo o la función de costo, es un ejemplo solamente, por lo cual la función puede ser modificada para una situación real de optimización. Pero cabe mencionar que la definición de ella es la parte más importante para aplicar estos algoritmos.

#Creamos una función a optimizar
#Función Objetivo o Función de Costo

costo.tiempo<-function(r){
 precio.total=0
 ultimo.arrivo=0
 salida.temprana=24*60
 
 for(i in 1:(length(r)/2)){
 #nombre=gente[[1]][i]
 origen=gente[[2]][i]
 Destino=c('LGA')
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 precio.total=precio.total+salida[,5]+regreso[,5]
 
 if(ultimo.arrivo<minutos(salida[,4])){ultimo.arrivo=minutos(salida[,4])}
 if(salida.temprana>minutos(regreso[,3])){salida.temprana=minutos(regreso[,3])}
 
 }
 tiempo_espera=0 
 for( i in 1:(length(r)/2)){
 origen=gente[[2]][i]
 Destino=c('LGA')
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 
 tiempo_espera=tiempo_espera+(ultimo.arrivo-minutos(salida[,4]))
 tiempo_espera=tiempo_espera+(minutos(regreso[,3])-salida.temprana)
 }
 if(ultimo.arrivo>salida.temprana){precio.total=precio.total+50} 
 
 return(precio.total+tiempo_espera)
 }

Si se hace un prueba con esta función se nota que devuelve un número entero, el cual representa el valor de la solución o configuración que asignamos a los viajantes.

#Ejemplo del uso de la función de costo
#Definimos una posible solución
 vector
 [1] 6 7 4 2 6 7 4 2 6 2 6 8
> costo.tiempo(vector)
[1] 3806

Paso 4

Para comparar métodos uso el método de búsqueda aleatoria. Consiste básicamente en tomar posibles soluciones y comparar entre las que se van encontrado, así con un proceso de comparación se elige la solución que minimiza  la función objetivo. Como es de esperar no es un método eficiente, pero permite acercarse a las ideas de la optimización estocástica.

#Búsqueda Aleatoria

optimizacion.aleatoria<-function(funcost){
 #dominio va ser igual al número máximo que puede tomar cada entrada
 best=999999999
 bestr=list()
 for( i in (0:1000)){
 r=sample(10,size=12,replace=TRUE) 
 cost=funcost(r)
 
 if(cost<best){
 best=cost
 bestr=r}
 }

 return(bestr)}

Cabe mencionar que el método fue adaptada para el problema, el algoritmo es estándar. Se ponemos a prueba el método y se evalúa la función de costo.

#Resultados del método Búsqueda Aleatoria

> optimizacion.aleatoria(costo.tiempo)
 [1] 8 9 6 7 2 8 8 7 6 8 7 8

> solucionba=c(8, 9, 6, 7, 2, 8, 8, 7, 6, 8, 7, 8)

> costo.tiempo(solucionba)
[1] 3689

> impri.horarios(solucionba)
[1] "Seymour" "BOS" "17:11-18:30 108" "18:25-20:34 205"
[1] "Franny" "DAL" "13:54-18:02 294" "15:07-17:21 129"
[1] "Zooey" "CAK" "8:27-10:45 139" "16:35-18:56 144"
[1] "Walt" "MIA" "17:07-20:04 291" "15:07-17:21 129"
[1] "Buddy" "ORD" "14:22-16:32 126" "16:35-18:56 144"
[1] "Les" "OMA" "15:03-16:42 135" "16:35-18:56 144"

Ahora usamos el algoritmo de recorrido simulados, tienen su origen en un problema de física. Por lo cual algunos de sus parámetros prevalecen en uso por analogía del fenómeno físico.

La idea es iniciar con una solución aleatoria y con un parámetro de temperatura alto, este último representa en la búsqueda la guía de un proceso de “enfriamiento” donde conforme decrece se van comparando los valores de la función objetivo. De los cuales se elige una posible solución al óptimo global que se espera obtener al tener un nivel de temperatura bajo.

El método es una adaptación al método Metropolis-Hastings, por lo cual en el algoritmo aparece una exponencial para comparar un valorar aleatorio y continuar en las comparaciones de los valores de la función objetivo.

#Algoritmo recorrido simulado
annealingoptimize<-function(costf,T=10000.0,cool=0.95,step=1){
 #Creamos una solución aleatória
 vec=sample(9,size=12,replace=TRUE)

 while(T>0.1){
 #Elejimos un índice
 i=sample(8,size=1)
 #Definimos la dirección del algortimos
 x=c(-1,0,1)
 dir=sample(x,1) 
 vecb=vec
 vecb[i]=vecb[i]+dir
 
 if(vecb[i]<0){vecb[i]=0}
 if(vecb[i]>9){vecb[i]=9}
 
 ea=costf(vec)
 eb=costf(vecb)
 
 
 if(eb<=ea){
 vec=vecb
 }
 else
 { 
 a=runif(1,min=0,max=1)
 
 if(a<exp((-eb-ea)/T))
 {
 vec=vecb
 
 } 
 else
 {
 vec=vec}
 
 
 }
 
 T=T*cool
 
 }
 return(vec)
 
 }

Hacemos una prueba de la optimización por este método. Los resultados fueron los siguientes.

#Método de Recorrido Simulado

annealingoptimize(costo.tiempo,T=10000.0,cool=0.95,step=1)
 [1] 8 7 4 7 8 3 7 3 1 5 8 7

> solucionrs=c(8,7,4,7,8,3,7,3,1,5,8,7)

#Evaluamos nuestra función de costo
> costo.tiempo(soucionrs)
[1] 4573
#Imprimimos nuestros vuelos

> impri.horarios(soucionrs)
[1] "Seymour" "BOS" "17:11-18:30 108" "15:07-17:21 129"
[1] "Franny" "DAL" "10:30-14:57 290" "15:07-17:21 129"
[1] "Zooey" "CAK" "17:08-19:08 262" "9:31-11:43 210" 
[1] "Walt" "MIA" "15:34-18:11 326" "9:31-11:43 210" 
[1] "Buddy" "ORD" "6:05-8:32 174" "12:31-14:02 234"
[1] "Les" "OMA" "16:51-19:09 147" "15:07-17:21 129"

Por último hago un pequeño código de un algoritmo genético, este fue adaptado para el problema que se está resolviendo. Lo que se debe de mencionar son las etapas que conforman dicho algoritmo: Mutación, Cruza, Generación de Población, Elitismo y procesamiento.

Cada etapa tienen varias opciones para ser realizada, la mutación puede ser solo una variación o afectación a la cadena de entrada. La cruza puede tener varios métodos para implementarse, y la idea base es “pegar” dos cadenas para obtener una nueva. La Población, tal cual como su nombre es, genera una cantidad de cadenas obtenidas de manera aleatoria y por “elitismo” elegir las mejores, es decir por un criterio se discrimina entre las cadenas.

Después de tener una Población inicial y haber elegido entre ellas,  se usa la mutación y la cruza para ir “mejorando” la población y por ultimo se  aproxima al óptimo.

El código que muestro es una versión simple y hasta cierto punto fea o mala, ya que hay pasos que pueden hacerse de manera más simplificada o compacta. Pero preferí mostrar todos los pasos por simples que fueran, pensando dar visibilidad en lugar de elegancia.

#Algoritmo Genético

genetic.optimize<-function(costf,popsize=50,step=1,mutprod=0.2,elite=0.2,maxiter=100){
 
 #Mutación
 mutacion<-function(vec){
 i=sample(12,size=1)
 a=runif(1,min=0,max=1)
 if((a<0.5)){
 if(i==1){return(c(vec[1]-step,vec[2:12]))}
 else if (i==12){return(c(vec[1:11],vec[12]-step))} 
 else
 {
 return(c(vec[1:(i-1)],vec[i]-step,vec[(i+1):12]))
 }
 }
 else if(vec[i]<9) 
 { 
 if(i==1){return(c(vec[1]+step,vec[2:12]))}
 else if (i==12){return(c(vec[1:11],vec[12]+step))}
 else 
 {
 return(c(vec[1:(i-1)],vec[i]+step,vec[(i+1):12]))
 }
 
 }
 
 }

 #Cruce entre dos cadenas de soluciones
 crossover<-function(r1,r2){
 i=sample(12,size=1) 
 if(i==1){return(c(r1[1],r2[2:12]))}
 else if (i==12){return(c(r1[1:11],r2[12]))}
 else 
 {return(c(r1[1:(i-1)],r2[i:12]))}
 }
 
 #Creación de la población inicial
 A=matrix(nrow=popsize,ncol=12)
 
 for(i in (1:popsize)){
 A[i,]=rbind(sample(9,size=12,replace=TRUE))
 }

 #Ganadores por generación
 topelite=as.integer(elite*popsize)
 
 #Loop principal
 for( i in (1:maxiter)){
 score=sapply(1:nrow(A),function(p){costf(A[p,])})
 B=cbind(A,score)
 Bdata<-data.frame(B)
 
 Bdata<-Bdata[with(Bdata,order(score)),]
 A=as.matrix(Bdata[1:topelite,1:12])
 while (length(A[,1])<popsize){
 a=runif(1,max=1,min=0)
 if(a<mutprod){
 c3=sample(topelite, size=1)
 A=rbind(A,mutacion(A[nrow=c3,1:12]))
 }
 else
 {
 c1=sample(topelite, size=1)
 c2=sample(topelite, size=1)
 A=rbind(A,crossover(A[nrow=c1,1:12],A[nrow=c2,1:12]))
 }
 } 
 
 }
 print(A[1,1:12])
 
 }

Se  compara el resultado y  se obtenemos lo siguiente.

#Método de Algoritmo Genético

 sol5
 6 5 4 5 6 7 4 5 5 7 5 5

costo.tiempo(sol5)
 2898

impri.horarios(sol5)
[1] "Seymour" "BOS" "13:40-15:37 138" "12:31-14:02 234"
[1] "Franny" "DAL" "10:30-14:57 290" "12:31-14:02 234"
[1] "Zooey" "CAK" "13:40-15:38 137" "15:07-17:21 129"
[1] "Walt" "MIA" "11:28-14:40 248" "12:31-14:02 234"
[1] "Buddy" "ORD" "12:44-14:17 134" "15:07-17:21 129"
[1] "Les" "OMA" "12:18-14:56 172" "12:31-14:02 234"

Es muy curioso como los resultados de los algoritmo muestran comportamientos muy distintos en cuando a las elecciones de vuelos y precios. Claro que esto es dependiente de la función de costo definida. El costo por medio de la solución de algoritmos genéticos es mucho mejor que las otras dos. Es casi 1000 unidades menos que la mejor de las otras dos.

Espero esta entrada de una idea de lo que es en general la optimización y el tipo de algoritmos que existen. Por su puesto que hay muchos temas que ni siquiera menciono, pero espero puedan consultarlos y que lo visto en esta entrada les ayude a tener una visión global del tema.

Comentario: El algoritmo que actualmente se emplea para hacer Machine Learning con grandes volúmenes de datos se llama Gradiente Descendiente Estocástico, se ha verificado que es adecuado para la mayoría de algoritmo de clasificación y predicción de ML, lo cual lo ha vuelto popular y se ha implementado en varias librerías, para ver  una breve, pero buena explicación el método recomiendo el vídeo de 12 minutos de Andrew Ng: https://class.coursera.org/ml-003/lecture/107

Referencias.

Libros:

http://www.amazon.com/Optimization-Machine-Learning-Information-Processing/dp/026201646X

http://www.princeton.edu/~sbubeck/Bubeck14.pdf

http://www.amazon.es/Stochastic-Optimization-Scientific-Computation-Schneider/dp/3540345590

Centros de investigación:

http://www.cs.manchester.ac.uk/mlo/

http://research.microsoft.com/en-us/groups/mlo/

Notas, artículos y cursos:

http://stoprog.org/SPTutorial/SPTutorial.php

http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc_print.pdf

http://online.stanford.edu/course/convex-optimization-winter-2014

http://www.econ.uiuc.edu/~roger/research/conopt/coptr.pdf