Machine learning desde cero (III)
Clasificación lineal y regresión logística

Repasito sobre la regresión

En los posts anteriores vimos cómo es posible predecir ciertos valores “de salida” a partir de ciertos valores “de entrada” a través de descubrir tendencias (rectas, planos, hiperplanos… con cierta pendiente) en observaciones anteriores.

A este cálculo de una posible tendencia lineal lo llamamos regresión lineal y vimos, de pasada, que es posible ajustar no sólo formas geométricas lineales sino también curvas (este tema lo trataremos en profundidad en futuros posts).

Desde el punto de vista del Machine Learning, puesto que la obtención o ajuste del patrón, modelo o tendencia se hace a partir de una familia de datos de entrenamiento en que los valores de salida son conocidos, la regresión lineal entra dentro de la familia de algoritmos de aprendizaje supervisado.

En este posts vamos a introducir otro de los algoritmos fundamentales de aprendizaje supervisado: el de la clasificación lineal.

La clasificación

εκστα ση… εκστα νο…
Chimo Bayo

En la regresión nuestro objetivo era predecir una variable respuesta en forma de de número real a partir de ciertas observaciones, por ejemplo, podíamos predecir precios de bienes inmuebles. En la clasificación, en cambio, nuestro objetivo es predecir un valor discreto, o dicho de otro modo, clasificar observaciones en categorías distintas: animal o vegetal, enfermo o sano, etc.

Las aplicaciones de la clasificación son muchas y muy útiles, entre las más conocidas están, por ejemplo, el filtro de spam (clasificar emails como spam o no spam) o la resolución de problemas de computer vision como el análisis de emociones (clasificar caras como triste, feliz…).

La clasificación lineal

Vamos a empezar por la clasificación lineal binaria, es decir, donde la variable respuesta que queremos predecir sólo puede tener dos categorías posibles (por ejemplo: spam o no spam, enfermo o sano, etc.). Existen diversos algoritmos para tratar de delimitar linealmente las categorías “de respuesta”, en este post usaremos el de la regresión logística.

Como veremos, la función objetivo a optimizar en la clasificación lineal como problema de regresión logística no tendrá una solución “en forma cerrada” (closed-form solution) o, dicho de otro modo, que no vamos a poder encontrar una fórmula con solución exacta.

Esta dificultad nos va a permitir introducir un algoritmo de búsqueda de soluciones de forma iterativa que, de hecho, aplicaremos a la resolución de muchos otros problemas de optimización que surgen en machine learning, incluyendo también el de la regresión lineal de los posts anteriores.

Empecemos por presentar una serie de datos de ejemplo en que se ha observado si una serie de pacientes presentaban o no presentaban una enfermedad cardíaca y la edad y frecuencia cardíaca máxima de cada paciente.

Nuestro objetivo es construir un modelo predictivo que nos permita predecir con cierto nivel de confianza si un nuevo paciente está sano o enfermo a partir de su edad y de medir su frecuencia cardíaca máxima. Hasta aquí vemos que el planteamiento es muy similar al de la regresión como predicción.

Una forma habitual de presentar gráficamente el problema de la clasificación lineal es a través de una nube de puntos sobre un plano en que los dos ejes cartesianos representarán nuestras variables predictivas (la Edad y la Frecuencia cardíaca máxima en nuestro ejemplo) y un color o forma distinta para cada punto que represente nuestra variable “objetivo” a predecir (en el ejemplo, si está Enfermo usaremos la cruz, si está sano la redonda). Veámoslo:

A ojo, podemos ver que los pacientes enfermos (las cruces) tienden a acumularse en el lado superior derecho, mientras que los sanos (los círculos) parecen acumularse en el lado inferior izquierdo. Da la sensación de que podríamos trazar una recta que separase estas dos áreas del gráfico. ¡Esta es precisamente la recta que vamos a aprender a calcular!

Pero Mitxel, ¿por qué no utilizar una curva en lugar de una recta? ¡Ciertamente! Es del todo posible utilizar curvas para clasificar, pero recuerda que este es un artículo sobre clasificación lineal, por lo que por ahora nos limitaremos a utilizar rectas (y planos e hiperplanos).

Pues bien, ¿qué modelo utilizar? Dado que nuestro objetivo es dividir el área en dos secciones, necesitaremos un modelo de recta que, al comparar la posición de los puntos respecto a ésta, arroje algún tipo de resultado numérico que informe de si el punto está a un lado u a otro.

Una primera idea podría ser utilizar los cuadrados de las distancias verticales de un punto a una recta, como hicimos en la regresión lineal, pero es fácil ver que el cuadrado (al perder su signo) no nos dice nada sobre si un punto cae a un lado u otro de una recta. Por otro lado, las distancias verticales (siguiendo el eje cartesiano vertical), tampoco parecen ser la mejor opción para determinar lo que está a un lado u otro de una recta (¿y si nuestra recta hipotética fuese completamente vertical?).

¡Pero la intuición era buena! En lugar de las distancias verticales, usaremos las distancias perpendiculares a nuestra recta hipotética, de este modo las (mínimas) distancias punto-recta y el signo de esta distancia nos informarán de si un punto está a un lado u otro. La distancia de cada uno de los puntos a la recta también será importante, puesto que en nuestro modelo un punto muy cercano a la recta, por definición, podría haber caído fácilmente a un lado u otro, mientras que los puntos más alejados serán los que clasificaremos con mayor confianza.

Así pues, para una recta, del tipo $ w_1x_1 + w_2x_2 + w_0 = 0 $ o, generalizando $w^Tx = 0$, donde $w_1$ y $w_2$ representan los weights o componentes horizontal y vertical a calcular, la distancia mínima o perpendicular a la recta será:

Por trigonometría básica y siendo $v$ la distancia vertical entre el punto y la recta:

Construyendo un triángulo rectángulo de componentes $w_1$ (horizontal, opuesto a $\alpha$) y $w_2$ (vertical, adyacente a $\alpha$):

Y combinando las dos idas:

Que es una fórmula muy bonita porque al haber desaparecido $w_2$, tiene forma generalizada (aplicable a dos o más dimensiones).

Ahora que tenemos una fórmula para las distancias, tenemos que encontrar una expresión optimizable, es decir, que nos permita encontrar los valores de los weights que mejor clasifican las observaciones realizadas en nuestra tabla (o cuan bien predicen futuras observaciones, que viene a ser lo mismo).

Una idea interesante podría ser, por ejemplo, utilizar la función signo para clasificar observaciones/generar predicciones:

Nótese que hemos prescindido del denominador de la expresión de la distancia $d$, puesto que al ser una constante positiva, no tiene ningún efecto sobre el resultado de la función signo, lo que nos permite simplificar.

Deberíamos intentar entonces maximizar que el número de observaciones $y_i$ que tengan el mismo signo que $\hat{y}_i$, o lo que es lo mismo, tratar de minimizar el número de casos del tipo:

La condición $y_i\,(w^T x_i) < 0$ sólo nos dice, para un punto, si está bien o mal clasificado. Para construir una función optimizable podríamos, por ejemplo, usar el sumatorio de la la función indicadora $\mathbf{1}[\cdot]$, que vale $1$ cuando la condición se cumple y $0$ cuando no.

El concepto de esta función cuya optimización representa una minimización de coste es popularmente conocida como loss function. Habitualmente no puede ser optimizada/resuelta de forma exacta, por lo que debemos encontrar métodos iterativos para resolverla.

Debido a que hemos elegido como base de nuestra loss function la función signo (que es analíticamente terrible, porque no es derivable), la única forma que tendríamos para tratar de optimizar nuestra función objetivo sería ploteando montones de rectas paralelas y con diferentes orientaciones en el espacio y calculando para cada una la cantidad de errores, hasta encontrar una (o varias) soluciones que resulten en el mínimo número de errores.

De forma visual:

Podemos escribir, de hecho, un script que minimice nuestra loss function de forma iterativa.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy as np


def minimo_errores_lineal(X, y):
    n = len(y)
    mejor_errores = n + 1
    mejor_frontera = None

    for k in range(3600):
        theta = np.pi * k / 3600.0
        normal = np.array([np.cos(theta), np.sin(theta)])
        proyecciones = X.dot(normal)
        orden = np.sort(proyecciones)
        cortes = np.concatenate([[orden[0] - 1],
                                 (orden[:-1] + orden[1:]) / 2.0,
                                 [orden[-1] + 1]])
        for t in cortes:
            prediccion = (proyecciones >= t).astype(int)
            errores = min((prediccion != y).sum(), ((1 - prediccion) != y).sum())
            if errores < mejor_errores:
                mejor_errores = errores
                mejor_frontera = (normal, t)
    return mejor_errores, mejor_frontera


EDAD = np.array([67, 65, 62, 63, 52, 56, 55, 61, 40, 53, 51, 41, 43, 57, 48,
                 64, 60, 61, 63, 37, 59, 56, 44, 52, 57, 54, 48, 66, 49, 58])
FC = np.array([108, 129, 160, 147, 155, 142, 148, 132, 154, 168, 160, 172, 160, 112, 139,
               144, 160, 173, 150, 187, 178, 153, 173, 162, 174, 160, 168, 158, 171, 162])

ENFERMO = np.array([0] * 15 + [1] * 15)

X = np.column_stack([EDAD, FC])
errores, frontera = minimo_errores_lineal(X, ENFERMO)

print("mejor frontera =", frontera)  # -> (array([0.8800633 , 0.47485639]), np.float64(120.86659558228276))
print("mínimo número de errores =", errores)  # -> 3

Así, una solución posible para nuestra frontera solución será:

Gráficamente:

Aunque nuestro primer modelo ingenuo nos ha permitido encontrar una solución para dos dimensiones, es muy importante notar que para el caso general nuestra solución es computacionalmente terrible, ya que al ir agregando variables/dimensiones, el número de orientaciones y paralelas que deberemos trazar crecerá exponencialmente:

Si echamos la vista atrás veremos que nuestros “problemas” comenzaron con la elección de una función con un muy pobre comportamiento analítico (la función signo). Una idea mucho más interesante será utilizar la función logística.

La función logística

La función logística tiene la siguiente forma estándar:

Su gráfica es una muy característica curva en forma de “S”, cuyo rango va de 0 a 1.

Usar esta función tiene diversas ventajas respecto a la función signo. En primer lugar, nos permitirá asignar a cada punto un valor decimal entre 0 y 1, lo que immediatamente podemos interpretar como una probabilidad. En segundo lugar, la curva es derivable en todo su dominio, y la derivada no es particularmente complicada.

Derivamos con la regla de la cadena (la derivada de $u^{-1}$ es $-u^{-2}u’$, y la derivada de $1+e^{-z}$ es $-e^{-z}$):

Ahora reescribimos la expresión para que aparezca la propia $\sigma$, separando el cuadrado del denominador en dos factores:

El primer factor es exactamente $\sigma(z)$. Y el segundo factor no es más que $1 - \sigma(z)$, como puede comprobarse:

Por tanto:

El modelo de regresión logística

La idea es realizar ahora el mismo trabajo que hicimos con la función signo, pero con la función logística.

Así pues, tomaremos del mismo modo la expresión “simplificada” de la distancia $w^T x$ y la “inyectaremos” en lugar de $z$.

Esta expresión, además, la podemos interpretar ahora como la probabilidad de que la observación pertenezca a la categoría “enfermo”.

Optimizar, ahora, será equivalente a encontrar los weights que determinan la recta que hace que nuestras observaciones sean las más probables posible.

La función objetivo: máxima verosimilitud

Para cada paciente, nuestro modelo asigna la probabilidad $p_i$ de estar enfermo y, por tanto, $1 - p_i$ de estar sano. Aprovechando que $y_i$ sólo puede tomar valores $0$ ó $1$, podemos expresar para cada punto la probabilidad de estar correctamente clasificado como un producto:

Nótese que si $y_i = 1$ queda $p_i$; si $y_i = 0$ queda $1 - p_i$.

De esta forma, la probabilidad de observar toda nuestra tabla de datos (aka verosimilitud) se puede expresar como el productorio de las probabilidades individuales:

Como maximizar un producto de muchos numeritos pequeños es numérica y algebraicamente incómodo, podemos aplicar el truco clásico de tomar los logaritmos. El logaritmo es una función creciente, de modo que el $w$ que maximiza $L(w)$ es el mismo que maximiza $\log L(w)$, y gracias a las propiedades de los logaritmos no sólo convertiremos el productorio en sumatorio, sino que también “bajaremos” los exponentes de las expresiones de probabilidad:

Como nuestro objetivo es, por convención, encontrar una loss function a minimizar, podemos invertir el signo de la expresión de máxima verosimilitud, lo que resulta en nuestra función objetivo definitiva, conocida como cross-entropy loss, logistic loss o log-loss:

Minimizar $J(w)$ equivale a maximizar la verosimilitud. Además, esta función es convexa, así que tiene un único mínimo global. ¡Muy buenas noticias!

Sin solución cerrada: el descenso de gradiente

En los posts anteriores derivábamos la función objetivo, igualábamos a cero y despejábamos $w$ con una fórmula cerrada. Si intentamos hacer lo mismo con $J(w)$, calcularemos su derivada e igualaremos a cero, pero nos toparemos con que la $\sigma$ (con su exponencial por dentro) hace imposible despejar $w$ de forma analítica. La solución es buscar el mínimo de forma iterativa (análogamente a como lidiamos anteriormente con la loss function con forma de función indicadora).

Aplicaremos aquí el famoso algoritmo para encontrar mínimos conocido como descenso de gradiente o gradient descent, que puede entenderse intuitivamente como a un excursionista perdido en la niebla que quiere bajar de la montaña: aunque no puede ver ningún valle a lo lejos, en cada punto puede notar hacia dónde baja más la pendiente y dar un pasito en esa dirección. Repitiendo el proceso muchas veces, acabará llegando a un valle (o, al menos, a uno de los valles).

Matemáticamente, “la dirección de máxima pendiente ascendente” es el gradiente (el vector de derivadas parciales), así que para bajar damos pasos en la dirección opuesta al gradiente:

donde $\eta$ (eta) es la tasa de aprendizaje (learning rate), un numerito que controla el tamaño del paso: si es demasiado grande, corremos el riesgo de pasarnos de largo y no converger nunca; si es demasiado pequeño, tardaremos una eternidad en llegar. Repetimos la actualización hasta que $w$ apenas cambie de una iteración a la siguiente.

El gradiente de la función objetivo

Sólo nos falta una pieza: la derivada de $J(w)$ respecto a $w$. Aquí es donde entra en juego la expresión de la derivada de la función logística (1) que calculamos anteriormente.

Primero, dos abreviaturas para no perdernos. Llamamos $z_i = w^T x_i$ al marcador lineal (proveniente de la idea de distancia) y $p_i = \sigma(z_i)$ a la probabilidad o sigmoide. Como $J$ depende de $w$ a través de $p_i$, que a su vez depende de $z_i$, aplicaremos la regla de la cadena:

donde $J_i = -\bigl[y_i \log p_i + (1 - y_i)\log(1 - p_i)\bigr]$ es el sumando $i$-ésimo.

En cuanto al primer término, derivamos respecto a $p_i$, recordando que la derivada de $\log(1-p)$ es $-\tfrac{1}{1-p}$:

El segundo término es la derivada de una sigmoide, cuya fórmula obtuvimos anteriormente en (1):

Nótese que si evaluamos el producto de estos dos términos podemos cancelar los denominadores:

\(\frac{\partial J_i}{\partial p_i}\cdot\frac{\partial p_i}{\partial z_i} = \frac{\partial J_i}{\partial z_i}\) \(= \left[-\frac{y_i}{p_i} + \frac{1 - y_i}{1 - p_i}\right] p_i(1 - p_i)\) \(= -y_i(1 - p_i) + (1 - y_i)\,p_i\)

Desarrollando y simplificando, casi todo se cancela:

En cuanto al tercer término, puesto que $z_i = w^T x_i$

En definitiva:

De donde se obtiene la fórmula “de libro”:

Que, además, puede expresarse en forma matricial, si $X$ es nuestra matriz de variables predictoras, $T$ el vector de valores objetivo y aplicamos la sigmoide elemento a elemento sobre $Xw$:

Nótese que, de hecho, es exactamente la misma forma que el gradiente de la regresión lineal $X^T(Xw - T)$, solo que ahora la predicción $Xw$ pasa por la sigmoide. La regla de actualización del descenso de gradiente queda entonces:

Manos a la obra: implementación en python

Nuestras variables predictoras viven en escalas muy distintas (la Edad ronda los 40–67 y la Frecuencia cardíaca los 108–187). Para realizar el descenso de gradiente es necesario normalizar las variables antes de entrenar: a cada una le restamos su media y la dividimos por su desviación típica. Una vez calculados los pesos en el espacio normalizado, los reconvertimos al espacio original para poder dibujar la frontera sobre nuestra nube de puntos.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import numpy as np


class LogisticRegression(object):
    def __init__(self, number_of_predictive_variables):
        self.k = number_of_predictive_variables
        self.w = np.zeros(number_of_predictive_variables + 1)
        self.mu = None
        self.sd = None

    @staticmethod
    def sigmoid(z):
        return 1.0 / (1.0 + np.exp(-z))

    def _design_matrix(self, x_data):
        x_normalized = (np.column_stack(x_data) - self.mu) / self.sd
        n = x_normalized.shape[0]
        return np.hstack([np.ones((n, 1)), x_normalized])

    def train(self, t_data, *x_data, learning_rate=0.1, iterations=20000):
        raw = np.column_stack(x_data).astype(float)
        self.mu = raw.mean(axis=0)
        self.sd = raw.std(axis=0)

        X = self._design_matrix(x_data)
        T = np.array(t_data, dtype=float)

        self.w = np.zeros(self.k + 1)
        for _ in range(iterations):
            gradient = X.T.dot(self.sigmoid(X.dot(self.w)) - T)
            self.w = self.w - learning_rate * gradient

    def predict_proba(self, *x_values):
        x = (np.array(x_values, dtype=float) - self.mu) / self.sd
        return self.sigmoid(np.concatenate([[1.0], x]).dot(self.w))

    def predict(self, *x_values):
        return int(self.predict_proba(*x_values) >= 0.5)

    def boundary_weights(self):
        # devuelve [b0, b1, b2] en el espacio ORIGINAL, tal que b0 + b1*x1 + b2*x2 = 0
        b0 = self.w[0] - np.sum(self.w[1:] * self.mu / self.sd)
        b_rest = self.w[1:] / self.sd
        return np.concatenate([[b0], b_rest])


AGE   = np.array([67, 65, 62, 63, 52, 56, 55, 61, 40, 53, 51, 41, 43, 57, 48,
                  64, 60, 61, 63, 37, 59, 56, 44, 52, 57, 54, 48, 66, 49, 58])
HRATE = np.array([108, 129, 160, 147, 155, 142, 148, 132, 154, 168, 160, 172, 160, 112, 139,
                  144, 160, 173, 150, 187, 178, 153, 173, 162, 174, 160, 168, 158, 171, 162])
SICK  = np.array([0]*15 + [1]*15)  # 15 sanos, 15 enfermos

model = LogisticRegression(2)
model.train(SICK, AGE, HRATE)

print("pesos (espacio normalizado) =", model.w)
print("pesos de la frontera (espacio original) =", model.boundary_weights())
print("prob. enfermo (edad 60, fc 175) =", model.predict_proba(60, 175))
print("prob. enfermo (edad 60, fc 120) =", model.predict_proba(60, 120))
# pesos (espacio normalizado) = [-0.30846899  1.79764042  3.37494467]
# pesos de la frontera (espacio original) = [-41.81646452   0.2237823    0.18840724]
# prob. enfermo (edad 60, fc 175) = 0.9898666714935286
# prob. enfermo (edad 60, fc 120) = 0.003077186040739993

Así, nuestra nueva solución usando este modelo será:

Nótese que cada iteración $T$ tiene un coste de \(O(n \cdot d)\), por lo que el coste del algoritmo será:

Puesto que $d$ aquí no es un exponente (a diferencia de antes), vemos que este nuevo modelo escalará mucho mejor con el número de dimensiones, una de sus ventajas principales.

Fíjate, además, en un detalle interesante: la frontera logística tiene un error más en el entrenamiento (4 vs 3) respecto a la frontera calculada a partir de la función indicatriz. Esto es porque no minimiza el conteo de errores, sino la entropía cruzada, que reparte la “culpa” de forma suave entre todos los puntos, en lugar de perseguir el mínimo exacto de errores.

Lo verdaderamente potente es que, al igual que ocurría con la regresión lineal, todo esto generaliza sin esfuerzo: con tres variables predictoras la frontera sería un plano, y con muchas más, un hiperplano que ya no podríamos dibujar pero que el descenso de gradiente calcularía exactamente igual y, además, en un tiempo razonable.

Próximamente

En los próximos posts exploraremos qué hacer cuando las categorías son más de dos (clasificación multiclase con la función softmax) y cómo el descenso de gradiente que acabamos de conocer se convierte en la herramienta que hace funcionar, ni más ni menos, que a las redes neuronales.