Vicente Rodríguez

Oct. 19, 2018

Logistic Regression desde cero con python

La regresión logistica es un modelo de machine learning que se usa para clasificar datos.

La regresión logistica es un modelo muy importante si quieres entender las redes neuronales y es lo primero que deberías aprender ya que tiene muchos conceptos importantes que se reutilizan en las redes neuronales.

Recomiendo leer el post sobre los arboles de decisión ya que explico muchos conceptos que se usan en este y muchos más modelos (como las redes neuronales).

Los datos

En este tutorial trabajaremos con datos de personas con diabetes. En total son 768 registros y 9 columnas, 8 de ellas son información sobre las personas (lo que llamamos features, en los datos tenemos caracteristicas como presión sanguinea, insulina, edad, numero de embarazos) y la ultima columna es la respuesta (1, 0) si la persona tiene diabetes o no (lo que llamamos target, en nuestros datos se llama outcome), lo que queremos es construir un modelo que nos diga si la persona tiene diabetes o no a partir de los datos (features).

Partes del modelo

El modelo cuenta con varias funciones, cada una tendra su sección para explicarlas mejor.

Función de probabilidad


Z = WX + b

La primera función a tomar en cuenta es la de probabilidad, esta función dará un resultado en porcentaje de que tan probable es que, dependiendo de los datos (X), la clase sea 1 o 0, en nuestro caso que tan probable es que la persona tenga diabetes o que no tenga diabetes. Como podemos notar esta función tiene tres partes importantes, W, X y b.

W

W (en Inglés weights), son los pesos que se le otorgaran a las caracteristicas(features: edad, presión sanguinea, insulina).

Los pesos asignan valores a cada característica dependiendo de la importancia de los datos a la hora de predecir si la persona tiene diabetes o no y como veremos más adelante cuando el modelo este en la fase de entrenamiento entendera cuales son las caracteristicas más importantes para tomar en cuenta cuando una persona tiene diabetes. Un ejemplo seria:


W = [[0.3, 0.1, 0.6, 0.8, 0.4, 0.6, 0.7, 0.4]]

Como podemos notar W es una matriz de una fila y 8 columnas, son 8 columnas porque tenemos 8 features.

X

X son los datos que el modelo usa para aprender, en este caso tenemos 8 columnas y 768 filas, tenemos que usar matrices para representar los datos y para realizar las operaciones usaremos multiplicación de matrices y aunque suene complicado es algo sencillo y cuando usemos la librería numpy solo necesitaremos una línea de código.

un ejemplo de x:


x = [

 [6, 148, 72,  35, 0, 33.6, 0.627, 50],

 [1, 85, 66, 29, 0, 26.6, 0.351, 31],

 [8, 183, 64, 0, 0, 23.3, 0.672, 32]

]

Aqui tenemos las 3 primeras filas con sus 8 columnas de los datos de personas con diabetes.

b

b (en Inglés bias), es un parametro que se añade a la función de probabilidad. Si graficaramos el modelo de regresión logistica lo que tendriamos es una línea dividiendo los datos en dos clases como en la siguiente imagen

logistic Regression

lo que el parametro b hace es mover la línea que divide los datos de su punto de origen, si no tuvieramos el parametro b la línea siempre comenzaría desde en 0, gracias a b podemos elegir desde que punto la linea comienza, en el ejemplo la línea comienza a proximadamente en 3, si este valor fuese mayor la línea comenzaría

más arriba y si fuese menor comenzaria más abajo.


b = 0

Multiplicación de matrices

Para que dos matrices se puedan multiplicar, el numero de columnas de la primera matriz tiene que ser igual al numero de filas de la segunda matriz.


x = [

 [6, 148, 72,  35, 0, 33.6, 0.627, 50],

 [1, 85, 66, 29, 0, 26.6, 0.351, 31],

 [8, 183, 64, 0, 0, 23.3, 0.672, 32]

]

La matriz X es 3x8

3 filas y 8 columnas


W = [[0.3, 0.1, 0.6, 0.8, 0.4, 0.6, 0.7, 0.4]]

La matriz W es 1x8

1 fila y 8 columnas

Pareciera que no podemos multiplicar estas dos matrices pero lo que haremos es transponer la matriz X, esto quiere decir que cambiaremos las columnas por las filas y la matriz pasara de ser 3x8 a 8x3.

La matriz X pasará de ser:


x = [

 [6, 148, 72,  35, 0, 33.6, 0.627, 50],

 [1, 85, 66, 29, 0, 26.6, 0.351, 31],

 [8, 183, 64, 0, 0, 23.3, 0.672, 32]

]

a ser:


X = [  

       [6,     1,     8     ],

       [148,   85,    183   ],

       [72,    66,    64    ],

       [35,    29,    0.    ],

       [0,     0,     0     ],

       [33.6,  26.6,  23.3  ],

       [0.627, 0.351, 0.672 ] ,

       [50,    31,    32    ]

]

8 filas y 3 columnas, con esta operación ya podemos multiplicar ambas matrices y la matriz resultante sera 3x1.

Para multiplicar una matriz se empiza multiplicando cada valor de la primera fila de la matriz A (W) por cada valor de la primera columna de la matriz B (X) y luego sumarlos.

W solo tiene una fila y es la fila con la cual multiplicaremos todas las columnas de la matriz B (X)


[0.3, 0.1, 0.6, 0.8, 0.4, 0.6, 0.7, 0.4]



[

 6, 

 148, 

 72, 

 35, 

 0, 

 33.6, 

 0.627, 

 50

]



(6 * 0.3) + (148 * 0.1) + (72 * 0.6) + (35 * 0.8) + (0 * 0.4) + (33.6 * 0.6) + (0.627 * 0.7) + (50 * 0.4)

el resultado es 128.3989 y este resultado se anota en la primera fila de nuestra nueva matriz (3x1)


[

  [128.3989],

  [],

  [],

]

ahora tenemos que realizar el mismo proceso para las columnas 2 y 3 de la matriz B(X), el resultado es el siguiente:


[

  [128.3989],

  [100.2057],

  [86.3504],

]

Como mencione los pesos W le otorgan importancía a cada feature de nuestros datos, entre mas grande sea el valor de W mayor importancía tiene el feature pero como podemos notar tenemos valores muy diferentes en cada columna, desde valores decimales hasta valores grandes como 148 lo cual provoca que el modelo piense que esa caracteristica(feature) tiene mas importancía, lo que tenemos que hacer es convertir todos los datos a un rango de -1 a 1, de esta manera los todos los datos tendran un valor mas parecido y mas razonable a la hora de realizar las multiplicaciones, esto lo haremos más adelante en el código.

Como resumen tenemos que la función:


Z = WX + b

regresa un valor generalmente decimal el cual es la probabilidad de que una serie de caracteristicas pertenezca a la clase 1 o 0 pero lo que nosotros queremos es que el modelo nos diga si la clase es 0 o 1 y es donde entra en juego la siguiente función.

Función Sigmoid

La función sigmoid es la siguiente:

```

1 / (1 + e ^ - z)

```

y en gráfica

función sigmoid

Lo que hace esta función es tomar el resultado de la función de probabilidad y regresar un valor mas acorde a lo que estamos buscando, si el resultado de la función de probabilidad fue mayor a 1, la función sigmoid regresara un valor muy cercado a 1, si el resultado fue menor de -1 la función sigmoid regresara un valor muy cercano a 0, nosotros queremos que el modelo nos regrese valores lo más cercanos posibles a la clase a la que pertenece, si estamos analizando a la primera persona y esta tiene diabetes, queremos que la función sigmoid regrese un valor de .99 y esto indicaría que el modelo trabaja correctamente, pero si el la función sigmoid regresa un valor como 0.5 o menor esto quiere decir que el modelo no aprendio muy bien.

Función de costo

Tenemos que tener una metrica para saber que tan bueno es nuestro modelo, de esto se encarga la función de costo:


cost = (- 1 / m) * sum(y * log(A) + (1 - y) * (log(1 - A)))

la función tiene varios parametros

Lo que el modelo busca es que los arrays y y A sean iguales. Cuando el modelo trabaja correctamente el resultado de la función de costo sera muy bajo, si el modelo tiene un desempeño pobre el resultado sera muy alto, el siguiente paso es encontrar la manera de que la función de costo sea baja y esto se logra encontrando mejores parametros W y b ya que estos son los que modifican la función de probabilidad y son los que otorgan la importancia a las caracteristicas(features).

Gradient Descent

Para encontrar los mejores parametros W y b tendremos que recurrir a calculo y derivadas, este es un tema que puede ser difícil de comprender si no tienes experiencía con calculo pero no es imposible, en resumen queremos encontrar los mejores valores para W y para b.

gif

Este gif es una buena explicación sobre lo que hace el algoritmo de gradient descent. Si graficamos la función de costo esta al ser una función convexa se dibuja con forma de medio ovalo en vertical y tiene un punto mínimo (En el gif esta representado por una estrella), lo que el algoritmo hace es buscar este punto mínimo ya que es donde la función se minimiza (este proceso también se conoce como optimización) y la forma de lograrlo es dando pequeños pasos, primero se empieza en un punto random (en el gif se empieza en la esquina superior izquierda) y conforme se van cambiando los valores de W y b el punto azul va trazando un camino que llega al punto mínimo.

Tengo un tutorial donde explico los detalles de este algoritmo y del backpropagation para las redes neuronales pero los conceptos también aplican para la regresión logistica. Ya en el tutorial sobre redes neuronales deberian de quedar claro estos conceptos.

Código

Crearemos el modelo desde cero para poder ver todos los pasos que lleva crear el modelo.

Recomiendo ver la libreta de Jupyter donde esta todo el código paso a paso y con ejemplos.

Link a la libreta

Empezaremos exportando las librerias necesarias:


from sklearn.model_selection import train_test_split

import numpy as np

import pandas as pd

Aquí esta el link para descargar los datos (si estas usando la libreta ya tiene los datos descargados).


features = dataset.drop(["Outcome"], axis=1)

X = np.array(features)

y = np.array(dataset["Outcome"])

Tenemos que tener las dos variables X, y donde estan las caracteristicas de las personas(X) y donde esta la clase a la que corresponde esa persona (y).

Un punto importante es normalizar los datos para que todos esten en un rango de -1 a 1 y no existan problemas.


X = (X - np.min(X)) / (np.max(X) - np.min(X))

Ahora pasaremos a crear funciones que nos ayudaran a crear el modelo

La función sigmoid:


def sigmoid(z):

   return 1 / (1 + np.exp(-z))

La función de probabilidad y el costo:


def forward_propagation(W, b, X, y):

  data_number = X.shape[0]

  Z = np.dot(W, X.T) + b

  A = sigmoid(Z)



  cost = (- 1 / data_number) * np.sum(y * np.log(A) + (1 - y) * (np.log(1 - A)))



  return A, cost

esta función regresa el resultado de la función de costo y el resultado de la función de probabilidad (A) mas adelante usaremos A para realizar el gradient descent, por ahora realizaremos el backpropagation.


def backward_propagation(X, A, y):

  data_number = X.shape[0]

  dW = (1 / data_number) * np.dot((A - y), X)

  db = (1 / data_number) * np.sum(A - y)



  return dW, db

((A - y)X es la derivada de W y (A - y) es la derivada de b)


def optimize(W, b, X, y, num_iterations, learning_rate):

  costs = []



  for i in range(num_iterations):

    A, cost = forward_propagation(W, b, X, y)

    dW, db = backward_propagation(X, A, y)



    W = W - learning_rate * dW

    b = b - learning_rate * db



    if i % 100 == 0:

      costs.append(cost)

      print ("Costo e iteracion %i: %f" % (i, cost))



  params = {

    "W": W,

    "b": b

  }



  gradients = {

    "dW": dW,

    "db": db

  }



  return params, gradients, costs

Esta función es mas grande que las anteriores porque aquí empezamos a juntar todo el código, hay partes importantes de esta función:


for i in range(num_iterations):

Usamos un ciclo for ya que queremos que el gradient descent se ejecute varias veces, esto porque se necesita llegar al punto mínimo, nosotros podemos indicar cuantas iteraciones realizara (num_iterations) y también cual es la taza de aprendizaje (learning_rate) si tenemos una taza de aprendizaje alta necesitaremos menos iteraciones y si la taza de aprendizaje es baja necesitaremos más iteraciones, estos dos parametros se modifican dependiendo del problema.


W = W - learning_rate * dW

b = b - learning_rate * db

en este paso se realiza el gradient descent, aquí el learning_rate multiplica al gradiente de W y de b para saber que tanto se movera el punto azul (como se puede ver en el gif de la sección anterior).


def predict(W, b, X):

  data_number = X.shape[0]

  y_prediction = np.zeros((1, data_number))



  Z = np.dot(W, X.T) + b

  A = sigmoid(Z)



  for i in range(A.shape[1]):

    y_prediction[0, i] = 1 if A[0, i] > 0.5 else 0



  return y_prediction

Una vez el modelo aprendio de los datos tendremos nuevos valores para W y para b, la función predict se encarga de usar estos nuevos valores para predecir las clases de nuevos datos


y_prediction[0, i] = 1 if A[0, i] > 0.5 else 0

en esta línea de código podemos ver como se asignan las clases dependiento de la probabilidad que obtuvieron, si la probabilidad es mayor a 0.5 la clase es 1, si es igual o menor la clase es 0

Recordemos que A*es el resultado de la función sigmoid y es un array parecido a:


[0.04, 0.99, 0.7, 0.5, 0.88, 0.33]

y cada elemento de este array se convertira en una clase ya sea 1 o 0.


def model(X_train, y_train, X_val, Y_val, num_iterations=2000, learning_rate=0.5):

  dimensions = X.shape[1]

  W = np.zeros(shape=(1, dimensions))

  b = 0



  params, gradients, costs = optimize(W, b, X_train, y_train, num_iterations, learning_rate)



  W = params["W"]

  b = params["b"]



  y_prediction_validation = predict(W, b, X_val)

  y_prediction_train = predict(W, b, X_train)



  print("train accuracy: {} %".format(100 - np.mean(np.abs(y_prediction_train - y_train)) * 100))

  print("test accuracy: {} %".format(100 - np.mean(np.abs(y_prediction_validation - y_val)) * 100))

para terminar el modelo la función model es la que encapsula todas las funciones anteriores y las ejecuta. También se encarga de imprimir que tan bueno es el modelo usando dos set de datos, el de entrenamiento y el de validación.


X_train, X_val, y_train, y_val = train_test_split(X, y, random_state=0, test_size=0.20)

con esta función dividimos los datos que tenemos en los sets de entrenamiento y validación.


model(X_train, y_train, X_val, y_val, num_iterations = 1000, learning_rate = 0.003)

Por ultimo podemos entrenar el modelo, pueden cambiar los parametros num_iterations y learning_rate para ver que tanto cambia, con estos parametros los resultados son los siguientes:


train accuracy: 64.00651465798046 %

test accuracy: 69.48051948051949 %