En esta libreta vamos a desarrollar los algoritmos de regresión logística, y vamos a aplicar los métodos a dos conjuntos de datos, uno donde se aplica directamente la regresión logística y otro donde se ejemplifica el uso de la regularización y los clasificadores polinomiales.
Igualmente vamos a probar el uso del descenso de gradiente, y el uso de herramientas de optimización que provienen de los paquetes de scipy.
In [1]:
%matplotlib inline
from scipy.stats import logistic
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image # Esto es para desplegar imágenes en la libreta
La función logística está dada por
$$ g(z) = \frac{1}{1 + e^{-z}}, $$la cual tambien se conoce como función probit (probabilidad de un bit). La cual es importante que podamos calcular en forma vectorial. Si bien el calculo es de una sola linea, el uso de estas funciones auxiliares facilitan la legibilidad del código.
In [2]:
def logistica(z):
"""
Calcula la función logística para cada elemento de z
@param z: un ndarray
@return: un ndarray de las mismas dimensiones que z
"""
# Introduce código aqui (una linea de código)
#---------------------------------------------------
return logistic.cdf(z)
#---------------------------------------------------
# prueba que efectivamente funciona la función implementada
# si el assert es falso regresa un error de aserción (el testunit de los pobres)
assert (np.abs(logistica(np.array([-1, 0, 1])) - np.array([ 0.26894142, 0.5, 0.73105858]))).sum() < 1e-6
Para probar la función vamos a graficar la función logística en el intervalo [-5, 5]
In [4]:
z = np.linspace(-5, 5, 100)
print z
with plt.xkcd():
plt.plot( z, logistica(z))
plt.title(u'Función logística', fontsize=20)
plt.xlabel(r'$z$', fontsize=20)
plt.ylabel(r'$\frac{1}{1 + \exp(-z)}$', fontsize=26)
Una vez establecida la función logística, vamos a implementar la función de costo sin regularizar para la regresión logística, la cual está dada por
$$ J(\theta) = -\frac{1}{T} \sum_{i=1}^T \left[ y^{(i)}\log(h_\theta(x^{(i)}) + (1 - y^{(i)})\log(1 - h_\theta(x^{(i)}))\right], $$donde
$$ h_\theta(x^{(i)}) = g(\theta^T x_e^{(i)}), $$las cuales fueron ecuaciones revisadas en clase.
In [5]:
def costo(theta, x, y):
"""
Calcula el costo de una theta dada para el conjunto dee entrenamiento dado por y y x
@param theta: un ndarray de dimensión (n + 1, 1)
@param x: un ndarray de dimensión (T, n + 1) donde la primer columna son puros unos
@param y: un ndarray de dimensión (T, 1) donde cada entrada es 1.0 o 0.0
@return: un flotante con el costo
"""
T = x.shape[0]
z = logistica(x.dot(theta))
#------------------------------------------------------------------------
# Agregua aqui tu código
return -np.sum( y*np.log(z) + (1 - y)*np.log(1 - z ))/T
#------------------------------------------------------------------------
# Otra vez el testunit del pobre (ya lo calcule yo, pero puedes hacerlo a mano para estar seguro)
theta = np.ones((2,1))
x = np.array([[1, 10],
[1, -5]])
y1 = np.array([[1],
[0]])
y2 = np.array([[0],
[1]])
y3 = np.array([[0],
[0]])
y4 = np.array([[1],
[1]])
assert abs(costo(theta, x, y1) - 0.01) < 1e-2
assert abs(costo(theta, x, y2) - 7.5) < 1e-2
assert abs(costo(theta, x, y3) - 5.5) < 1e-2
assert abs(costo(theta, x, y4) - 2.0) < 1e-2
De la misma manera, para poder implementar las funciones de aprendizaje, vamos a implementar el gradiente de la función de costo. El gradiente de la función de costo respecto a $\theta$ es (como lo vimos en clase) el siguiente:
$$ \frac{\partial J(\theta)}{\partial \theta_j} = -\frac{1}{T} \sum_{i=1}^T (y^{(i)} - h_\theta(x^{(i)})x_j^{(i)} $$y a partir de las ecuaciones individuales de puede obtener $\nabla J(\theta)$, la cual no la vamos a escribir en la libreta para que revisen en sus notas como se puede resolver este problema en forma matricial. Si bien para el descenso de gradiente podemos utilizar directamente el gradiente negado, al implementar métodos de optimización avanzados, necesitamos que el gradiente sea efectivamente el gradiente.
In [6]:
def gradiente(theta, x, y):
"""
Calcula el gradiente del costo de la regrasión logística, para una theta, conociendo un conjunto de aprendizaje.
@param theta: un ndarray de dimensión (n + 1, 1)
@param x: un ndarray de dimensión (T, n + 1) donde la primer columna son puros unos
@param y: un ndarray de dimensión (T, 1) donde cada entrada es 1.0 o 0.0
@return: un ndarray de mismas dimensiones que theta
"""
T = x.shape[0]
#------------------------------------------------------------------------
# Agregua aqui tu código
G = logistica(x.dot(theta))
result = (1./T)*(x.T.dot(G - y))
return result
#------------------------------------------------------------------------
# Otra vez el testunit del pobre (ya lo calcule yo, pero puedes hacerlo a mano para estar seguro)
theta = np.ones((2, 1))
x = np.array([[1, 10],
[1, -5]])
y1 = np.array([[1],
[0]])
y2 = np.array([[0],
[1]])
y3 = np.array([[0],
[0]])
y4 = np.array([[1],
[1]])
assert abs(0.00898475 - gradiente(theta, x, y1)[0]) < 1e-4
assert abs(7.45495097 - gradiente(theta, x, y2)[1]) < 1e-4
assert abs(4.95495097 - gradiente(theta, x, y3)[1]) < 1e-4
assert abs(-0.49101525 - gradiente(theta, x, y4)[0]) < 1e-4
Ahora vamos a desarrollar las funciones necesarias para realizar el entrenamiento y encontrar la mejor $\theta$ de acuero a la función de costos y un conjunto de datos de aprendizaje.
Para este problema, vamos a utilizar una base de datos sintética proveniente del curso de Andrew Ng que se encuentra en coursera. Supongamos que pertenecemos al departamente de servicios escolares de la UNISON y vamos a modificar el procedimiento de admisión, y en lugar de utilizar un solo exámen (EXCOBA) y la información del cardex de la preparatoria, hemos decidido aplicar dos exámenes, uno sicométrico y otro de habilidades estudiantiles. Dichos exámenes se han aplicado el último año aunque no fueron utilizados como criterio. Así, tenemos un historial entre estudiantes aceptados y resultados de los dos exámenes. El objetivo es hacer un método de regresión que nos permita hacer la admisión a la UNISON tomando en cuenta únicamente los dos exámenes y simplificar el proceso. Recuerda que esto no es verdad, es solo un ejercicio.
Bien, los datos se encuentran en el archivo admision.txt el cual se encuentra en formato cvs (osea los valores de las columnas separados por comas. Vamos a utilizar el mismo método para cargar datos que el que utilizamos para regresión lineal, y vamos a graficar la información para entender un poco los datos.
In [11]:
datos = np.loadtxt('admision.txt', comments='%', delimiter=',')
x, y = datos[:,0:-1], datos[:,-1:]
x = np.c_[np.ones((x.shape[0], 1)), x]
plt.plot(x[y.ravel() == 1, 1], x[y.ravel() == 1, 2], 'sr', label='aceptados')
plt.plot(x[y.ravel() == 0, 1], x[y.ravel() == 0, 2], 'ob', label='rechazados')
plt.title(u'Ejemplo sintético para regresión logística')
plt.xlabel(u'Calificación del primer examen')
plt.ylabel(u'Calificación del segundo examen')
plt.axis([20, 100, 20, 100])
plt.legend(loc=0)
Out[11]:
Vistos los datos un clasificador lineal podría ser una buena solución.
Ahora vamos a implementar el método de descenso de gradiente, casi de la misma manera que lo implementamos para regresión lineal (por lotes)
In [12]:
def descenso_rl_lotes(x, y, alpha, epsilon=1e-4, max_iter=int(1e4), costos=False):
"""
Descenso de gradiente por lotes para resolver el problema de regresión logística con un conjunto de aprendizaje
@param x: un ndarray de dimensión (T, n + 1) donde la primer columna son puros unos
@param y: un ndarray de dimensión (T, 1) donde cada entrada es 1.0 o 0.0
@param alpha: Un flotante (típicamente pequeño) con la tasa de aprendizaje
@param epsilon: Un flotante pequeño como criterio de paro. Por default 1e-4
@param max_iter: Máximo numero de iteraciones. Por default 1e4
@param costos: Un booleano para saber si calculamos el historial de costos o no
@return: theta, costo_hist donde costo es ndarray de dimensión (n + 1, 1) y costo_hist es un
ndarray de dimensión (max_iter,) con el costo en cada iteración si costos == True, si no
regresa None
"""
T, n = x.shape[0], x.shape[1] - 1
theta = np.zeros((n + 1, 1))
costo_hist = np.zeros(max_iter) if costos else None
if costos:
for iter in xrange(max_iter):
theta += -alpha*gradiente(theta,x, y)
if np.linalg.norm(theta) <= epsilon:
return theta, costo_hist
costo_hist[iter] = costo(theta, x, y)
else:
for iter in xrange(max_iter):
theta += -alpha*gradiente(theta,x, y)
if np.linalg.norm(theta) <= epsilon:
return theta, costo_hist
#--------------------------------------------------------------
# Agregar aqui tu código
#
# Recuerda utilizar las funciones que ya has desarrollado
#--------------------------------------------------------------
return theta, costo_hist
Para probar la función de aprendizaje, vamos a aplicarla a nuestro problema de admisión. Primero recuerda que tienes que hacer una exploración para encontrar el mejor valor de alpha. Así que utiliza el código de abajo para ajustar alpha
In [18]:
alpha = 1e-4
mi = 500
_, costo_hist = descenso_rl_lotes(x, y, alpha, epsilon=1e-4, max_iter=mi, costos=True)
plt.plot(np.arange(mi), costo_hist)
plt.title(r'Evolucion del costo en las primeras iteraciones con $\alpha$ = ' + str(alpha))
plt.xlabel('iteraciones')
plt.ylabel('costo')
Out[18]:
Una vez encontrada la mejor $\alpha$, entonces podemos calcular $\theta$ (esto va a tardar bastante), recuerda que el costo final debe de ser lo más cercano a 0 posible, así que agrega cuantas iteraciones sean necesarias:
In [19]:
theta_lotes, _ = descenso_rl_lotes(x, y, alpha, max_iter = int(1e6))
print theta_lotes
print costo(theta_lotes, x, y)
Es interesante ver como el descenso de gradiente no es eficiente en este tipo de problemas, a pesar de ser problemas de optimización convexos.
Bueno, este método nos devuelve $\theta$, pero esto no es suficiente para ecir que tenemos un clasificador, ya que un método de clasificación se compone de dos métodos, uno para aprender y otro para predecir. Recuerda que para realizar la predicción con el método de regrasión lineal y el criterio MAP, no es necesario calcular la regrasión logística completa (lo que hace un método de predicción más rápido).
In [11]:
def predictor(theta, x):
"""
Predice los valores de y_hat (que solo pueden ser 0 o 1), utilizando el criterio MAP.
@param theta: un ndarray de dimensión (n + 1, 1)
@param x: un ndarray de dimensión (T, n + 1) donde la primer columna son puros unos
@return: y_hat un ndarray de dimensión (T, 1) donde cada entrada es 1.0 o 0.0
"""
#-------------------------------------------------------------------------------------
# Agrega aqui tu código sin utilizar la función logística
return np.where(x.dot(theta) > 0, 1, 0)
#--------------------------------------------------------------------------------------
¿Que tan bueno es este clasificador? ¿Es que implementamos bien el método?
Vamos a contestar esto por partes. Primero, vamos a graficar los mismos datos pero vamos a agregar la superficie de separación, la cual en este caso sabemos que es una linea recta. Como sabemos el criterio para decidir si un punto pertenece a la clase 1 o cero es si el valor de $\theta^T x_e^{(i)} \ge 0$, por lo que la frontera entre la región donde se escoge una clase de otra se encuentra en:
$$ 0 = \theta_0 + \theta_1 x_1 + \theta_2 x_2, $$y despejando:
$$ x_2 = -\frac{\theta_0}{\theta_2} -\frac{\theta_1}{\theta_2}x_1 $$son los pares $(x_1, x_2)$ de valores en la forntera. Al ser estos (en este caso) una linea recta solo necesitamos dos para graficar la superficie de separación.
In [12]:
x1_frontera = np.array([20, 100]) #Los valores mínimo y máximo que tenemos en la gráfica de puntos
x2_frontera = -(theta_lotes[0] / theta_lotes[2]) - (theta_lotes[1] / theta_lotes[2]) * x1_frontera
plt.plot(x[y.ravel() == 1, 1], x[y.ravel() == 1, 2], 'sr', label='aceptados')
plt.plot(x[y.ravel() == 0, 1], x[y.ravel() == 0, 2], 'ob', label='rechazados')
plt.plot(x1_frontera, x2_frontera, 'm')
plt.title(u'Ejemplo sintético para regresión logística')
plt.xlabel(u'Calificación del primer examen')
plt.ylabel(u'Calificación del segundo examen')
plt.axis([20, 100, 20, 100])
plt.legend(loc=0)
Out[12]:
Y para que tengas una idea de lo que debería de salir, anexo una figura obtenida con el código que yo hice:
In [13]:
Image(filename='ejemplo_logistica_1.png')
Bueno, ya vimos que el método del descenso de gradiente, si bien es correcto y fácil de implementar, es bastante ineficiente en cuanto a tiempo (inclusive para un problema bastante simple, por lo que vamos a utilizar la función minimizeque provee scipy la que vamos a utilizar con el algoritmo BFGS (explicado en clase en forma somera). Ejecuta la celda de abajo para revisar la documentación de minimize
In [14]:
from scipy.optimize import minimize
minimize?
Como puedes ver, el método BFGS es el método por default. Los parámetros que hay que agregar son:
x0funjacfun y jac en forma de tuplaoptions como un diccionario. Las opciones pueden ser:La función regresa un objeto resultado res, del cual solo nos interesa res.x, el resultado de theta.
Como $\theta$ en este método debe de ser un ndim de una sola dimensión, y nosotros hemos estado usando $\theta$ como un vector columna pues podemos modificar un poco estas funciones con unas nuevas (esto de hecho es posible hacerlo modificando las funciones originales, pero podría causar mas problemas de conceptos hacerlo). De hecho una vez funcionando, un buen científica de la computación (o un buen desarrollador) aplicaría una etapa de refactoring, esto es, modificar todo lo que obscurece y hace menos ineficiente, o menos estético el código.
Veras que si lo corres varias veces, en muchas ocasiones vamos a tener que el resultado no se pudo lograr por falta de precisión en la función del gradiente. Esto es debido a la forma en que estamos calculando el gradiente, y a que si la solución inicial realiza una partición pésima del estado, entonces hay problemas para estimar el inverso del hessiano, por problemas de estabilidad numérica del método. Por esta razón se inicializa el algoritmo con un valor de $\theta$ obtenido en forma aleatoria (valores pequeños).
Lo interesante es que el algoritmo avisa cuando no hay convergencia, y cuando hay se obtiene en muy pocas iteraciones (19 en promedio)
In [15]:
theta0 = 1e-2 * np.random.rand(x.shape[1])
print theta0
funcion = lambda theta, x, y: costo(theta.reshape(-1,1), x, y)
jacobiano = lambda theta, x, y: gradiente(theta.reshape(-1,1), x, y).ravel()
res = minimize(x0=theta0,
fun=funcion,
jac=jacobiano,
args = (x, y),
method='BFGS',
options= {'maxiter': 400, 'disp': True})
print res
Otra forma de obtener el resltado es utilizando un método tipo simplex que no necesita calcular ni el jacobiano, ni la inversión del hessiano. Que si bien requiere sensiblemente de más iteraciones, para un problema relativamente simple como este es suficiente (sin embargo para problemas de regresión lineal con muchos atributos, o redes neuronales, ya no será suficiente).
In [16]:
theta0 = np.zeros(x.shape[1])
print theta0
funcion = lambda theta, x, y: costo(theta.reshape(-1,1), x, y)
jacobiano = lambda theta, x, y: gradiente(theta.reshape(-1,1), x, y).ravel()
res = minimize(x0=theta0,
fun=funcion,
args = (x, y),
method='Nelder-Mead',
options= {'maxiter': 400, 'disp': True})
print res
Y ahora veamos la partición del espacio con este método, la cual es ligeramente diferente a la obtenida con el descenso de gradiente:
In [17]:
theta_NM = res.x.reshape(-1,1)
x1_frontera = np.array([20, 100]) #Los valores mínimo y máximo que tenemos en la gráfica de puntos
x2_frontera = -(theta_NM[0] / theta_NM[2]) - (theta_NM[1] / theta_NM[2]) * x1_frontera
plt.plot(x[y.ravel() == 1, 1], x[y.ravel() == 1, 2], 'sr', label='aceptados')
plt.plot(x[y.ravel() == 0, 1], x[y.ravel() == 0, 2], 'ob', label='rechazados')
plt.plot(x1_frontera, x2_frontera, 'm')
plt.title(u'Ejemplo sintético para regresión logística')
plt.xlabel(u'Calificación del primer examen')
plt.ylabel(u'Calificación del segundo examen')
plt.axis([20, 100, 20, 100])
plt.legend(loc=0)
Out[17]:
In [ ]: