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 clasificació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 clasificación utilizando regresión logística. A continuación en la Sección 3, se compararán diferentes clasificadores utilizando sklearn para resolver un problema de clasificación real.
El modelo de regresión linear discutido en la práctica anterior se propone bajo la asunción de que la variable respuesta es cuantitativa. En muchas situaciones, la variable respuesta es cualitativa o categóricas, por ejemplo, tipos de manzanas: golden, reineta, Royal Gala. El proceso de asociar cada manzana a un tipo concreto se conoce como clasificación. Por otro lado, habitualmente los métodos utilizados para clasificación estiman, en primer lugar, la probabilidad de pertenencia a cadad una de las posibles categorías, y utilizan estas probabilidades para decidir la clasificación final. Este procedimiento es parecido a la aproximación utilizada por los métodos de regresión lineal.
La regresión logística no trata de predecir el valor numérico de una variable cuantitativa, sino que intenta predecir la probabilidad de que una observación dada pertenezca a una cierta clase (categoría). Asumiendo que estamos, por razones de simplicidad y sin perdida de generalidad, en un caso binario únicamente con dos clases, la regresión logística da como resultado la siguiente probabilidad: $$ P\left[ \ y = 1\ |\ \mathbf{x}\ \right] $$
de tal forma que la probabilidad de pertenencia a la otra clase es: $$ P\left[\ y = 0\ |\ \mathbf{x}\ \right] = 1 - P\left[\ y = 1\ |\ \mathbf{x}\ \right] $$
Por lo que la salida de nuestro modelo de regresión logística debe estar comprendido en el lintervalo $[0,1]$.
La regresión logística asume que podemos utilizar una frontera lineal para separar el espacio de observaciones en dos regiones, una por cada clase. Supongámos que tenemos observaciones en un espacio bi-dimensional, esto es, las observaciones poseen dos características cada una $\mathbf{x} = [x_1 x_2]^T$. En la figura siguiente se pueden apreciar observaciones pertenecientes a dos clases, una representada con círculos negros y otra con círculos blancos.
En la misma figura se presentan tres fronteras lineales diferentes, $H_1, H_2$ y $H_3$. Para el caso de que las observaciones posean más dimensiones, hablamos de hyperplanos frontera, como se puede observar en la siguiente figura
Las fronteras lineales, en dos dimensiones, están definidas por la ecuación de una recta, por lo que puede ser expresada como:
$$h(\mathbf{x}) = w_0 + w_1 x_1 + w_2 x_2$$que expresado en notación vectorial sería:
$$\mathbf{w}^T \mathbf{x} = 0$$(Nota: el alumno avisado puede relacionar rápidamente esta ecuación con el modelo de regresión lineal presentado en la práctica previa.)
Suponga a hora el alumno que se tiene una nueva observación $\mathbf{x}^* = [x^*_1,x_2^*]$. Cuando evaluamos nuestra frontera lineal con esta nueva observación se pueden dar tres posibilidades:
$\mathbf{x}^*$ se encuentra en la región de la clase positiva ($y = 1$). De esta forma $h(\mathbf{x}^*) > 0$. Cuanto mayor sea el valor del $|h(\mathbf{x}^*)|$, más lejos se encontrará el punto de la frontera.
$\mathbf{x}^*$ se encuentra en la región de la clase negativa ($y = 0$). De esta forma $h(\mathbf{x}^*) < 0$. Cuanto mayor sea el valor del $|h(\mathbf{x}^*)|$, más lejos se encontrará el punto de la frontera.
$\mathbf{x}^*$ se encuentra en exactamente en la frontera. De esta forma $h(\mathbf{x}^*) = 0$.
Sin embargo, tal y como hemos dicho, necesitamos convertir este valor de $h(\mathbf{x})$ en una probabilidad, es decir en un valor entre $[0,1]$, que además conserve las propiedades de $h(\mathbf{x})$. En concreto, si estamos estimando la probabilidad de pertenecer a la clase positiva ($y = 1$), que el valor sea más cercano a 1 cuanto mayor sea el valor de $h(\mathbf{x})$, que sea $0.5$ cuando $h(\mathbf{x}) = 0$ y que sea cercano a 0 cuanto menor sea $h(\mathbf{x})$. El objetivo actual es, pues, econtrar una transformación $\theta(z)$, tal que:
$$P\left[ \ y = 1\ |\ \mathbf{x}\ \right] = \theta\left(h(\mathbf{x})\right) = \theta\left(\mathbf{w}^T\mathbf{x}\right)$$y que cumpla con todos los requisitos que hemos impuesto anteriormente. Existen varias funciones $\theta(z)$ que servirían a nuestro propósito. Una de ellas, la empeleada para la regresión logística es la que se conoce como función sigmoide o función logística, cuya expresión es:
$$\theta(z) = \frac{1}{1 + e^{-z}}$$Ya debería estar un poco más claro por qué se le llama regresión logística a este modelo.
En el siguiente código veremos la representación de la función logística.
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def logistic(z):
return 1.0 / (1.0 + np.exp(-z))
z = np.arange(-7, 7, 0.1)
phi_z = logistic(z)
plt.plot(z, phi_z)
plt.axvline(0.0, color='k')
plt.ylim(-0.1, 1.1)
plt.xlabel('z')
plt.ylabel('$\phi (z)$')
# y axis ticks and gridline
plt.yticks([0.0, 0.5, 1.0])
ax = plt.gca()
ax.yaxis.grid(True)
plt.tight_layout()
plt.show()
Recordando que $z = \mathbf{w}^T\mathbf{x}$, podemos modificar la forma de la regresión logística cambiando los pesos $\mathbf{w}$, que recordamos son los que vamos a estimar minimizando alguna función de coste.
De esta forma ya tenemos un modelo de regresión logística de forma que la probabilidad de que una observación $\mathbf{x}$ pertenezca a la clase 1 viene dada por:
$$P\left[ \ y = 1\ |\ \mathbf{x}\ \right] = \frac{1}{1+e^{-\mathbf{w}^T\mathbf{x}}}$$y, consecuentemente, diremos que $\mathbf{x} \in C_1$ sii $P\left[ \ y = 1\ |\ \mathbf{x}\ \right] > 0.5$.
Una vez que hemos definido nuesto modelo, falta saber cómo estimamos los coeficientes del modelo $\mathbf{w}$. En regresión lineal, utilizando el modelo de regresión lineal ($\mathbf{w}^T\mathbf{x}$) más un ruido Gaussiano de media cero y varianza $\sigma^2$, llegamos a una expresión cerrada para la estimación de los coeficientes del modelo:
$$\hat{\mathbf{w}} = \left(X^TX\right)^{-1}X^T\mathbf{y}$$Por contra, en el caso de regresión logística no existe una expresión cerrada, vamos a tener que utilizar otros métodos. En la sección siguiente definimos la función objetivo que optimizaremos, en el caso de regresión logística será el likelihood, y describiremos uno de los métodos más utilizados para su maximización, el gradiente ascent.
El método que utilizaremos para estimar los parámetros del modelo de regresión logística será el de maximizar la función de verosimilitud:
$$\hat{\mathbf{w}} = \underset{\mathbf{w}}{\operatorname{argmax}} L(\mathbf{x},\mathbf{w}) $$donde $L(\mathbf{x},\mathbf{w})$ es la función de verosimilitud del modelo de regresión logística. Por lo tanto, el trabajo ahora reside en definir la función de verosimilitud. En este punto podemos pensar que nuestro modelo de regresión logística, para una obervación dada $x_n$, definido por: $P\left[ \ y_n = 1\ |\ \mathbf{x_n}\ \right] = \theta\left(h(\mathbf{x_n})\right) = \theta\left(\mathbf{w}^T\mathbf{x_n}\right)$, no es nada más que la probabilidad de que y = 1. Podemos por lo tanto pensar en un experimento en el que lanzamos una modena, y asignamos una probabilidad $p$ a que salga cara ($y = 1$), y por lo tanto, una probabilidad $1-p$ a que salga cruz ($y = 0$). La función de masa de probabilidad que define el lanzamiento de una modena es:
$$P(y|x) = p^y\left(1-p\right)^{(1-y)}$$Este modelo probabilístico es completamente análogo al nuestro donde $p = P(y_n = 1|\mathbf{w},\mathbf{x_n})$ y $1-p = P(y_n = 0|\mathbf{w},\mathbf{x_n})$, de esta forma, la función de verosimilitud para una observación dada $(\mathbf{x}_n,\mathbf{y}_n)$ sera:
$$L(\mathbf{x_n},\mathbf{w})= P(y_n|\mathbf{w},\mathbf{x_n}) = P(y_n = 1|\mathbf{w},\mathbf{x_n})^{y_n}P(y_n = 0|\mathbf{w},\mathbf{x_n})^{(1-y_n)}$$Donde:
$$P(y_n = 0|\mathbf{w},\mathbf{x_n}) = 1 - P(y_n = 1|\mathbf{w},\mathbf{x_n}) = 1 - \frac{1}{1+e^{-\mathbf{w}^T\mathbf{x}_n}} = \frac{e^{-\mathbf{w}^T\mathbf{x}_n}}{1+e^{-\mathbf{w}^T\mathbf{x}_n}}$$Asumiendo que las observaciones se toman de forma independiente, podemos escribir la función de verosimilitud conjunta, para todas las observaciones como:
$$L\left(\mathbf{X},\mathbf{w}\right) = p(\mathbf{y}|\mathbf{w},\mathbf{X}) = \prod_{n = 1}^{N} P(y_n = 1|\mathbf{w},\mathbf{x_n})^{y_n}P(y_n = 0|\mathbf{w},\mathbf{x_n})^{(1-y_n)}$$Podemo escribir de forma un poco más compacta esta expresión utilizando la siguiente nomenclatura $P(y_n = 1|\mathbf{w},\mathbf{x_n}) = \rho_n$, donde $\rho_n$ es la probabilidad de que la observación $n-ésima$ tenga como valor de la variables respuesta $y_n = 1$.
$$L\left(\mathbf{X},\mathbf{w}\right) = \prod_{n = 1}^{N} \rho_n^{y_n}(1-\rho_n)^{(1-y_n)}$$Durante el curso hemos visto que en estos casos parece razonable utilizar la función logaritmo de $L$ ($\log L$), lo cual simplifica los calculos posteriores; de esta forma:
$$\log\left(L\left(\mathbf{X},\mathbf{w}\right)\right) = \sum_{n = 1}^{N} y_n\log(\rho_n) + (1-y_n)\log(1-\rho_n)$$Nota: cabe destacar que la función objetivo $J(\mathbf{w}) = - \log\left(L\left(\mathbf{X},\mathbf{w}\right)\right)$ se conoce como Cross-Entropy Error, y que maximizar $\log\left(L\left(\mathbf{X},\mathbf{w}\right)\right)$ equivale a minimizar $J(\mathbf{w})$.
So far, so good. Ahora estamos en disposición de una función de verosimilitud (en realidad del logaritmo), que bien utilizada nos permitirá estimar los parámetros de mi modelo, esto es: vamos a utilizar como estimadores $\hat{\mathbf{w}}$, aquellos que maximizan la función $\log(L)$. Sin embargo hemos comenzado este lab indicando que no vamos a poder encontrar una expresión cerrada para $\hat{\mathbf{w}}$ como resultado de derivar la función $\log(L)$. Cómo podemos resolver este problema? Gradient Descent (Ascent) al rescate.
El descenso por gradiente es un método de optimización, que nos permite minimizar (o maximizar) funciones objetivo de forma iterativa; como puede ser $\log(L)$ o cross-entropy error. La premisa fundamental es que dado un punto, en este caso $\mathbf{w}$, del espacio de busqueda, es decir del dominio de la función optimizar, vamos actualizar (update) ese punto, en el siguiente paso, en la dirección opuesta al gradiente de la función a optimizar. Esto sería para minimizar, si queremos maximizar nos moveremos en la dirección del gradiente. Vamos a generar intuición analizando un caso muy sencillo, en el que queremos minimizar una función objetivo $J(w)$, donde el parámetro $w$ es un escalar. Si somos afortunados, la función $J(w)$ será convexa.
Vamos a ilustrar el método del gradient descent en esta función.
El primer paso es seleccionar un punto, $w^{(0)}$, en el espacio de búsqueda (dominio de la función $J(w)$).
A continuación, calculamos el gradiente de la función objetivo para minimizar en el punto $w^{(0)}$: $\nabla J\left(w^{(0)}\right)$. El grandiente tiene la virtud de indicar la dirección de mayor crecimiento de la función $J(w)$, de forma que si el objetivo es "viajar" en dirección al mínimo deberemos dirigirnos en la dirección "menos el gradiente".
De forma genérica, si estamos en el paso $\tau$, los pesos en el paso $\tau + 1$ se actualizan de acuerdo con la regla:
$$w^{(\tau + 1)} = w^{(\tau)} - \eta \nabla J\left(w^{(\tau)}\right) $$
$\eta$ nos indica el cuanto variamos $w$ en cada iteración. Si $\eta$ es muy pequeño, en cada iteración la variación de $w$ será muy pequeña, por lo que el algoritmo puede tardar mucho en converger, y en ocasiones puede no llegar nunca al mínimo. Mientras que si $\eta$ es muy grande podemos tener problemas, de forma que en vez de ir hacia el mínimo nos movamos en dirección contraria.
En la siguiente figura vemos un ejemplo en el que el espacio de búsqueda es bidimensional, $\mathbf{w} = [w_0,w_1]$
Bien munidos, como estamos, con una función de verosimilitud que maximizar y un procedimiento para maximizar funciones, sólo puedrá atribuirse a nuestra perece si no somos capaces de obtener un estimador para los modelos de nuestro parámetro. Para completar todos nuestros preparativos necesitamos calcular el gradiente de nuestra función objetivo. El alumno puede demostrar que:
$$\frac{\partial L(\mathbf{w})}{\partial w_i} = \sum_{n = 1}^{N} (y_n - \rho_n)x_{n_i} $$Si escribimos esta expresión en forma matricial:
$$\nabla L(\mathbf{w}) = X^T(Y - T)$$NOTA esta ecuación no incluye el término de bias, para este término, la expresión de la parcial es:
$$\frac{\nabla L(\mathbf{w})}{\partial w_0} = \sum_{n = 1}^{N}(y_n-\rho_n)x_{n_{0}} = \sum_{n = 1}^{N}(y_n-\rho_n)$$En esta práctica vamos a estudiar el problema de intentar predecir el llamado Apgar index. Este índice se utiliza en la práctica clínica para evaluar el estado general del recién nacido después del parto. En este índice es habitual hacer la siguiente división, si Apgar > 6 se considera una situación normal, si Apgar $\le$ 6 situación anormal o sospechosa. El problema que se nos plantea es el siguiente: somos capaces de predecir el índice Apgar utilizando diferentes parámetros de la señal de fetal heart rate (FHR). El FHR es una señal que se adquiere de forma rutinaria durante el trabajo de parto, por lo que sería muy útil si fuesemos capaces de, simplemente analizando esa señal, predecir cuál va a ser el índice Apgar, esto permitiría tomar medidas rápidas, o incluso adelantar el parto.
Nosotros vamos a abordar este problema desde el punto de vista de clasificación, de forma que los recien nacidos pueden pertenecer a una de las dos siguientes clases: Normal, $y = 1$, si Apgar > 7; y anormal $y = 0$, si Apgar $\le$ 7.
Las variables explicativas, características o features serán los parámetros estimados a partir del FHR, en concreto serán:
La fuente de esta base de datos es Dr. Diogo Ayres-de-Campos, Fac. Med. Univ. Porto
La primera parte de esta práctica consistirá en leer los datos utilizando Pandas y generar nuestra variable respuesta binaria, $y$.
In [40]:
#read data using pandas
import pandas as pd
import numpy as np
fhr = pd.read_csv("FHR-Apgar.csv", delimiter=';')
fhr
Out[40]:
In [56]:
# Hay que valorar si merece la pena eliminar toda una fila que contenga algún NaN o directamente todas las
# que contengan todos los NaN, aunque puede ser que nos interesen algunos valores de las columnas
fhr = fhr.dropna()
#Esta instruccion genera un array donde cualquier valor de Apgar 1 > 7, lo ponemos a 1, en otro caso a 0
fhr['Y'] = np.where(fhr['Apgar 1'] > 7, 1, 0)
fhr
Out[56]:
In [57]:
%matplotlib inline
x_var_name = list(fhr)[4:-1] #Matriz X de caracteristicas
y_var_name = list(fhr)[-1] #Vector Y de respuesta
# X --> | |
# | ML |
# Y --> | |
X = fhr[x_var_name]
Y = fhr[y_var_name]
Una vez tenemos claros nuestros datos el siguiente paso, siempre, siempre, debe ser separar nuestros datos en dos conjunto: train y test.
Eventualmente, el conjunto train, resultante después del split anterior, puede ser subdividido para la selección de hyperparameters, pero la división incial nos va a servir para evaluar las prestaciones de nuestro modelo final (una vez seleccionados los hyperparameters). Test set
Este tipo de particiones se pueden realizar sencillamente con sklearn Cross-validation: evaluation estimator performance
In [93]:
from sklearn.model_selection import train_test_split
#test_size es el porcentaje de datos que utilizamos para test y otros para datos
#random_state son numeros pseudo aleatorios para que se haga el split en distintos datos
#Esto devuelve 2 matrices X y 2 vectores Y
#Si modificamos el test_size demasiado (muy bajo), estamos siendo muy optimistas para nuestro modelo
#El objetivo es ajustar (20-30% de test) lo mejor posible nuestro modelo
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
#(156, 8)
#(67, 8)
Pareciera que estamos listos para aplicar nuestros algoritmos de clasificación .... Uff, todavía falta un pequeño detalle, que en ocasiones marca la diferencia. Tenemos la siguiente pregunta:
Si la respuesta a la anterior pregunta es Sí, entonces se recomienda normalizar (estandarizar) nuestros datos Normalize). Se puede probar la diferencia utilizando normalización y sin utilizarla.
Esta fase se puede realizar de forma sencilla con sklearn: Preprocessing data
In [92]:
from sklearn.linear_model import LogisticRegression
#cremos un modelo vacio
log_reg = LogisticRegression() #Habra que modificar sus parametros
#Ajustamos el modelo
log_reg.fit(X_train, y_train)
print log_reg
#Hacemos la prediccion
y_pred = log_reg.predict(X_test)
#Vemos la accuracy
from sklearn.metrics import accuracy_score
acc = accuracy_score(y_test, y_pred)
print "Soy el puto amo de ML, con un % de: ", acc, '%'
In [ ]:
Por fin llegó. Ahora estamos en condiciones de entrenar nuestro primer algoritmo de clasificación. Comenzaremos por uno de los más sencillos: Regresión Logística. Este modelo nos servirá para aprender la API que utiliza sklearn para todos los algoritmos de aprendizaje máquina Logistic Regression. La diferencia entre ellos estriba en el hecho de los parámetros que tiene cada algoritmo. Algunos dirán que en la selección de esos hyperparameters radica una de las partes artísticas del machine learning. A mi me gusta más hablar de la parte artesana, creo que somos más artesanos que artistas.
En lo que sigue: aprenderemos a ajustar un modelo de clasificación, que será el de regresión logística y a evaluar las prestaciones. A continuación compararemos con diferentes modelos, en los que tendremos que seleccionar diferentes hyperparemeters.
Manos a la obra
In [ ]:
Vamos a probar con otro tipo de clasificadores. Por ejemplo: naïve bayes, [Support Vector Machines] (http://scikit-learn.org/stable/modules/svm.html#svm-classification), Random Forest
In [97]:
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
#grid search, todas las posibles combinaciones y hacemos cross-validation, rejilla logaritmica desde 10⁻⁴ hasta 10¹.
# con 25 valores para cada uno
#Es otro tipo de modelado
tuned_parameters = {'gamma' : np.logspace(-4, 1, 25), 'C' : np.logspace(1, 4, 25)}
svm = SVC(kernel = 'rbf')
grid = GridSearchCV(svm, param_grid = tuned_parameters, cv = 10, n_jobs = -1)
grid.fit(X_train, y_train)
y_pred_svm = grid.predict(X_test)
acc2 = accuracy_score(y_test, y_pred_svm)
print "SVM algorythm: ", acc2, '%'