Métodos de Regresión

Authors

Óscar Barquero Pérez (oscar.barquero@urjc.es), Carlos Figuera Pozuelo (carlos.figuera@urjc.es) y Rebeca Goya Esteban (rebeca.goyaesteban@urjc.es)

27 de marzo de 2016


Este obra está bajo una licencia de Creative Commons Reconocimiento-NoComercial-CompartirIgual 4.0 Internacional.

1 Introducción a la práctica

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.

2 Regresión usando Regresión Lineal

2.1 Regresión Lineal

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.

2.2 Regresión Lineal Polinómica

Regresión Lineal con polinomios

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.

3. Comparación de métodos de regresión con sklearn

3.1 Base de datos de regresión

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:

  1. CRIM: per capita crime rate by town
  2. ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
  3. INDUS: proportion of non-retail business acres per town
  4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
  5. NOX: nitric oxides concentration (parts per 10 million)
  6. RM: average number of rooms per dwelling
  7. AGE: proportion of owner-occupied units built prior to 1940
  8. DIS: weighted distances to five Boston employment centres
  9. RAD: index of accessibility to radial highways
  10. TAX: full-value property-tax rate per $10,000
  11. PTRATIO: pupil-teacher ratio by town
  12. B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
  13. LSTAT: % lower status of the population

Ejercicio 1. Lectura de la base de datos

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()


CRIM       0
ZN         0
INDUS      0
CHAS       0
NOX        0
RM         0
AGE        0
DIS        0
RAD        0
TAX        0
PTRATIO    0
B          0
LSTAT      0
PRICE      0
dtype: int64
Out[38]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE
count 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000
mean 3.593761 11.363636 11.136779 0.069170 0.554695 6.284634 68.574901 3.795043 9.549407 408.237154 18.455534 356.674032 12.653063 22.532806
std 8.596783 23.322453 6.860353 0.253994 0.115878 0.702617 28.148861 2.105710 8.707259 168.537116 2.164946 91.294864 7.141062 9.197104
min 0.006320 0.000000 0.460000 0.000000 0.385000 3.561000 2.900000 1.129600 1.000000 187.000000 12.600000 0.320000 1.730000 5.000000
25% 0.082045 0.000000 5.190000 0.000000 0.449000 5.885500 45.025000 2.100175 4.000000 279.000000 17.400000 375.377500 6.950000 17.025000
50% 0.256510 0.000000 9.690000 0.000000 0.538000 6.208500 77.500000 3.207450 5.000000 330.000000 19.050000 391.440000 11.360000 21.200000
75% 3.647422 12.500000 18.100000 0.000000 0.624000 6.623500 94.075000 5.188425 24.000000 666.000000 20.200000 396.225000 16.955000 25.000000
max 88.976200 100.000000 27.740000 1.000000 0.871000 8.780000 100.000000 12.126500 24.000000 711.000000 22.000000 396.900000 37.970000 50.000000

Ejercicio 2. Separar nuestros datos en test y training

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


(354, 13)
(152, 13)

 Ejercicio 3

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 "------------------"


bias
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-74-8bbd048d18d6> in <module>()
     13 
     14 for n,c in zip(labels, coefs):
---> 15     print n,c(round(c, 3))
     16     print "------------------"

TypeError: 'numpy.float64' object is not callable

Ejercicio 4

En este punto vamos a realizar la evaluación del modelo. Para ello, se va a realizar la predicción de los valores de las casas sobre el conjunto de test y utilizar una métrica conveniente para regresión.


In [75]:
from sklearn.metrics import r2_score

y_pred = lin_reg.predict(X_test)

#r2 score

print r2_score(y_test, y_pred)


 0.673528086535

Ejercicio 5

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)


Usando SVM:  0.0965880270632
Usando Random Forest:  0.81041284888