Este obra está bajo una licencia de Creative Commons Reconocimiento-NoComercial-CompartirIgual 4.0 Internacional.
El objetivo de esta lab es que el alumno se familiarice con los conceptos básicos de python, en concreto de la utilización del módulo Sklearn para resolver problemas de regresión.
Este lab está organizado de la siguiente forma: en la Sección 2, se realiza una introducción más formal al proceso de regresión utilizando regresión lineal. A continuación en la Sección 3, se compararán diferentes modelos de regresión utilizando sklearn para resolver un problema de clasificación real.
Uno de los modelos de aprendizaje estadístico más simples, aunque no por ello menos potente, es el modelo de Regresión Lineal.
Como se ha visto en clase de teoría, cuando utilizamos modelos de regresión lineal se asume que el valor de la señal que queremos estimas (predecir), que se conoce como variable respuestas u objetivo, es una combinación lineal de las variables de entrada, también conocidas como variables explicativas. En este caso, el objetivo planteado como un problema de estimación es obtener el valor de los coeficientes que multiplican a las variables explicativas, $w_i$. En notación matemática:
$$\hat{y}= w_0 + w_1 x_1 + ... + w_p x_p = \mathbf{w^{T}x}$$Nomenclatura: los coeficientes $w_i$ se conocen también como pesos; es común que al coeficiencte $w_0$ se le llame bias.
En regresión lineal se eligen los coeficientes $\mathbf{w}$ de forma que se minimice el error cuadrático, es decir que el error entre el valor real de $y$ y el proporcionado por nuestro modelo $\hat{y} = \mathbf{w^{T}x}$ sea el menor posible, para todos los valores. Es decir buscamos resolver el siguiente problema de minimización:
$$\underset{w}{min\,} {\left|\left| X w - y\right|\right|_2^2}$$Este problema de minimización se puede resolver recurriendo al, nunca suficientemente bien ponderado, Principio de Ortogonalidad.
El objetivo es, pues, estimar el vector $\mathbf{y}$ con el vector que mejor lo aproxime y que esté en el subespacio definido por las columnas de $X$, que llamaremos $\mathbf{\hat{y} = w^{T}x}$. Esto lo conseguimos proyectando el vector $\mathbf{y}$ en dicho subespacio. Esto es equivalente a que el vector error, $\mathbf{ \varepsilon = y - \hat{y}}$ sea ortogonal al subespacio de $X$
Gracias a la profesora Inmaculada Mora Jiménez por la figura
De esta forma, si el error es perpendicular al subespacio de $X$, sabemos que:
$$X^{t}\varepsilon = 0$$Entonces
$$X^{T}\left(\mathbf{y} - X\mathbf{w}\right) \Rightarrow X^{T}\mathbf{y} = X^{T}X\mathbf{w}$$Despejando los pesos:
$$\mathbf{w} = \left(X^{T}X\right)^{-1}X^{T}\mathbf{y}$$Vemos que el Principio de Ortogonalidad nos ofrece una solución cerrada para la estimación de los pesos de nuestro modelo de regresión lineal. Culpa nuestra será, por pereza o impericia, que no seamos capaces de programar tan sencilla solución utilizando numpy.
Documentación: http://scikit-learn.org/stable/modules/linear_model.html#polynomial-regression-extending-linear-models-with-basis-functions
El modelo de regresión lineal, tal y como se ha planteado, sólo serviría para resolver problemas en los que la relación entre la variable respuesta y las variables explicativas es lineal, tanto en los coeficientes como en las variables. Cuando existe alguna relación no lineal entre $y$ y las variables explicativas se pude extender sencillamente nuestro modelo de regresión lineal.
Esta extensión es común a muchos algoritmos de aprendizaje estadístico. Esta extensión se consigue construyendo polinomios a partir de las variables explicativas (características). Así, por ejemplo, en un caso con dos variables explicativas, el modelo de regresión lineal básico sería:
$$\hat{y} = w_0 + w_1 x_1 + w_2 x_2$$Si el resutado no es un plano, sino que es un paraboloide, se puede plantear la siguiente extensión a polinomios de segundo orden:
$$\hat{y} = w_0 + w_1 x_1 + w_2 x_2 + w_3 x_1 x_2 + w_4 x_1^2 + w_5 x_2^2$$El alumno que observa por primera vez esta extensión suele pensar que ya no tratamos con un modelo lineal, pero en este punto se debe hacer notar que la linealidad está en los coeficientes, y el modelo se convierte en completamente lineal por el simple expediente de un cambio de variable:
$$z = [x_1, x_2, x_1 x_2, x_1^2, x_2^2]$$Este cambio de variable nos permite reescribir el modelo de regresión lineal como:
$$\hat{y} = w_0 + w_1 z_1 + w_2 z_2 + w_3 z_3 + w_4 z_4 + w_5 z_5$$Como se puede observar este modelo de regresión es de la misma forma que los que hemos resuelto hasta ahora, y por lo tanto se puede resolver con las mismas técnicas. Considerando potencias más elevadas se pueden obtener modelos con una flexibilidad mucho mayor, es decir, se pueden resolver problemas en los que la relación entre la variable respuesta y las explicativas es muy compleja. Sin embargo, se corre un grave riesgo aumentando la complejidad del modelo (en este caso la complejidad es equivalente al mayor orden de polinomio), puesto que aumentamos la probabilidad de que el modelo sobreajuste. Más sobre esto en próximas secciones.
A continuación se propone ajustar un modelo de regresión lineal con bases no lineales.
En esta práctica vamos a utilizar la base de datos de los precios de casas en Bostón. Podéis encontrar toda la información en el repositorio de Bases de Datos para Machine Learning UCI-Boston
En la siguiente lista se proporciona la información de las características que se tiene por cada casa.
Attribute Information:
En este ejercicio se pide que el alumno lea la base de datos utilizando pandas. En aras de la brevedad, solo vamos a realizar la lectura del fichero csv y crear un describe del DataFrame. Se recomienda encarecidamente a los alumnos que realicen una exploración completa de la base de datos.
In [72]:
#read data using pandas
import pandas as pd
import numpy as np
data = pd.read_csv("boston.csv")
In [38]:
#describe
print np.sum(data.isnull()) #Usando isnull nos proporciona un booleano en caso de que tenga algun NaN, por lo que
# si algun numero es distinto de 0 hay algun NaN
data.describe()
Out[38]:
Como siempre que estemos trabajando con modelos de machine learning. El primer paso, después de haber realizado el análisis explroratorio es separar en nuestros conjunto de test y training. Cabe recordar en este punto que esta separación es obligatoria para poder reportar una estimación de la performance de nuestro modelo lo más realista posible. En el caso de disponer de un menor número de datos, se puede optar por evaluar las prestaciones de los modelos mediante un esquema de cross-validation, pero en general este último esquema produce resultados un poco más optimistas. Este tipo de particiones se pueden realizar sencillamente con sklearn Cross-validation: evaluation estimator performance
Lo primero será obtener nuestra matrix de características $X$ y el vector respuesta $y$
In [73]:
from sklearn.model_selection import train_test_split
X = data[list(data)[:-1]]
Y = data[list(data)[-1]]
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.3, random_state=0)
print X_train.shape
print X_test.shape
En este ejercicio vamos entrenar nuestro primer modelo de regresión. En este caso, igual que para clasificación, os recomiendo por comenzar con el modelo más sencillo posible: Regresión Lineal. El modelo lineal permite aumentar nuestro conocimiento sobre la naturaleza del problema
In [74]:
from sklearn import linear_model
from sklearn.metrics import accuracy_score
lin_reg = linear_model.LinearRegression()
lin_reg.fit(X_train, y_train)
coefs = [lin_reg.intercept_]
coefs.extend(list(lin_reg.coef_))
labels = ["bias"]
labels.extend(list(data)[:-1])
for n,c in zip(labels, coefs):
print n,c(round(c, 3))
print "------------------"
In [75]:
from sklearn.metrics import r2_score
y_pred = lin_reg.predict(X_test)
#r2 score
print r2_score(y_test, y_pred)
En este punto, evolucionamos nuestro modelo y podemos utilizar modelos más complicados. Por ejemplo modelos no lineales: Vamos a compara con SVM y Random Forest
Igual que hicimos en el lab sobre clasificación, vamos a tener que utilizar una estrategia de búsqueda de los hiperparámetros de los modelos. Para ello vamos a utilizar GridSearch y cross-validation
In [82]:
# SVM
from sklearn.svm import SVR
from sklearn.model_selection import GridSearchCV
#Se hace GridSearch para modelos que tengan hiperparametros
svm = SVR()
busqueda = {'C':np.logspace(1, 4, 10), 'epsilon': np.logspace(-4, 1, 10)}
gridSearch = GridSearchCV(svm, param_grid = busqueda, cv = 10, n_jobs = -1)
gridSearch.fit(X_train, y_train)
y_pred_svr = gridSearch.predict(X_test)
print "Usando SVM: ", r2_score(y_test, y_pred_svr)
#RandomForest
from sklearn.ensemble import RandomForestRegressor
rand_forest = RandomForestRegressor()
busquedaRF = {'n_estimators':[20, 50, 500]}
gridSearch = GridSearchCV(rand_forest, param_grid = busquedaRF, cv = 10, n_jobs = -1)
gridSearch.fit(X_train, y_train)
y_pred_rf = gridSearch.predict(X_test)
print "Usando Random Forest: ", r2_score(y_test, y_pred_rf)