Regresión lineal y descenso de gradiente con Python

En machine learning, el objetivo principal es encontrar un modelo que explique el comportamiento de un sistema (en el amplio sentido de la palabra). A partir de unos datos de entrenamiento, un sistema de aprendizaje automático ha de ser capaz de inferir un modelo capaz de explicar, al menos en su mayoría, los efectos observados. Pero también aplicar ese aprendizaje. Por ejemplo, un sistema de machine learning muy lucrativo para las empresas anunciantes es aquél que dado un perfil de usuario (datos de entrada A), sea capaz de predecir si pinchará o no (salida B) sobre un anuncio publicitario de, por ejemplo, comida para gatos. No es sencillo crear un modelo capaz de predecir el comportamiento del usuario (o sí), pero en todo caso, existen diferentes técnicas que nos permiten abordar el problema. En el caso del ejemplo que acabamos de ver, el modelo debería ser capaz de clasificar a los usuarios en dos clases diferentes, los que pulsarán y los que no pulsarán el anuncio de comida de gatos. Otro tipo de problema es aquél en el que se pretende crear un modelo capaz de predecir un valor o una cantidad. Por ejemplo, si conozco los metros cuadrados, el número de habitaciones y el número de baños de una vivienda en venta, ¿seré capaz de predecir por qué precio se venderá? Uno podría hacerse una idea observando el precio de venta de alguna casa igual o parecida que se haya vendido antes. Si no tenemos ese dato, podremos extrapolar el precio comparando el precio de venta de otras casas y comparando esa casa con la nuestra. Intuitivamente sabemos que a más metros, más precio ¿pero cuanto más? Es probable que entre los metros de la casa y el precio haya una relación lineal, así que aprovechándonos de ello, es relativamente sencillo construir un modelo (lineal, claro) que nos ayude a predecir el precio de una casa. Así que manos a la obra.
Vamos a basarnos en una técnica estadística llamada regresión lineal, que conceptualmente es bastante sencilla, pero es muy usada en este tipo de problemas. Antes de seguir vamos a generar un conjunto de datos ficticios con los que empezar a trabajar.

import random
import matplotlib.pyplot as plt

def gen_data(n, bias, varianza):
    x = []
    y = []
    for i in range(0, n):
        x.append(i)
        y.append((i + bias) + random.uniform(0, 1) * varianza)
    return x, y

x,y = gen_data(100, 25, 50)
plt.scatter(x, y)

Este código genera una serie de datos que sigue básicamente una recta alrededor de la cual se crean los puntos de forma aleatoria usando una función de distribución uniforme con una varianza dada (50 en este caso). Echemos un vistazo a como quedan los datos.

Por seguir con nuestro ejemplo de las casas puede imaginar que en el eje de abscisas se representan los metros cuadrados que tiene la casa y en el las ordenadas el precio en miles de euros por el que se vendió la casa.
La idea es, a partir de todos estos datos, crear un modelo que sea capaz de predecir el precio de cualquier casa dados sus metros cuadrados. En este caso es muy obvio que un buen predictor es una línea recta (por eso lo de lineal) que se ajuste lo más posible a la nube de puntos. Así pues, nuestro problema ahora es encontrar esa recta, la de mejor ajuste.
La ecuación de una recta en el plano es la siguiente.
\[y=a+bx\]
Donde a se corresponde con el punto de corte de la recta con el eje de ordenadas y b es la pendiente de la recta, es decir, lo rápido que crece o decrece. La variable independiente x representaría en este caso los metros cuadrados de la casa y la variable y el precio de venta. Por lo tanto, si somos capaces de buscar unos buenos valores para a y b, tendremos un modelo lineal para predecir el precio.
Podemos hacerlo de diferentes formas. Típicamente, la recta de regresión se ha calculado a partir de métodos estadísticos, pero (siempre hay un pero) si el número de casas que tenemos es muy grande (no lo llaman big data por gusto) o en vez de una sola variable independiente tenemos varias (metros, número de habitaciones, número de baños, metros de jardín, etc.) puede llegar a ser tedioso cuando no inabordable. Así pues vamos a usar un método para ajustar la recta llamado descenso de gradiente (gradient descent en la lengua de Shakespeare). Pero vayamos por partes. Lo primero que necesitamos es una forma de medir lo bien o mal que se ajusta una recta a la nube de puntos.
A cada par de valores que probemos para a y b lo vamos a llamar hipótesis, y seguidamente comprobaremos cómo se ajusta esa hipótesis calculando un factor de error, que es la suma de lo que, de media, se equivoca dicha hipótesis (recta) al predecir cada uno de los puntos de la nube. Este error se calcula con.
\[e=\frac{1}{2m}\sum_{i=1}^{m}(h(x^{i})-y^{i})^{2}\]
Donde \(h(x^{i})\) es la hipótesis que queremos probar. Por ejemplo, si asignamos los valores \(a=1\) y \(b=2\), tenemos que \(h(x^{i})=1+2x^{i}\) siendo \(x^{i}\) cada uno de los valores de x para los metros cuadrados de las casas e \(y^{i}\) los precios de dicha casa. Como se puede observar, lo que estamos haciendo es restar al precio de la casa que predice nuestra hipótesis (recta) el valor de venta real. Esa diferencia es el error, que se va sumando para todas las casas de las que disponemos de datos. El resultado de la resta (error) se eleva al cuadrado para evitar que unos errores se compensen con otros, de forma que el error siempre sea positivo. Resumiendo, lo que buscamos son los valores de a y b (nuestra hipótesis) que minimice el error devuelto por esta función. Ahí es donde entra en juego nuestro algoritmo del descenso de gradiente.
El algoritmo de descenso de gradiente hace uso de las derivadas (sí, eso que en el instituto pensaste que jamás te serviría para nada) y más concretamente de las derivadas parciales (sí, eso que en la universidad pensaste que jamás te serviría para nada). La función de error que acabamos de ver es una función convexa (no voy a demostrarlo aquí), lo que es muy interesante, ya que sólo tiene un valor mínimo que es global a toda la función. Puedes imaginarla con una forma similar a esta.

Imagina que el punto A es nuestra hipótesis inicial, y como queremos encontrar el valor de a y b que minimicen la función de error, lo que se busca es llegar al punto B. Conceptualmente, la derivada nos indica la velocidad a la que crece o decrece una función en uno de sus puntos. Por lo tanto, lo que debemos hacer es ir moviendo el punto A en la dirección de mayor pendiente para llegar a B. Como hay dos variables en juego, a y b, debemos calcular las derivadas parciales de la función de error con respecto a las dos variables.
\[a=a-\alpha \frac{\partial }{\partial a}e(a,b)\]
\[b=b-\alpha \frac{\partial }{\partial b}e(a,b)\]
El parámetro \(\alpha\) es lo que se llama learning rate, y permite controlar la velocidad de descenso hacia el punto mínimo de la función. Un valor bajo hará que la convergencia sea lenta, mientras que un valor alto puede hacer que nos pasemos de frenada y acabemos dando santos de un lado a otro de la función sin alcanzar nunca el mínimo. Valores típicos son 0.1, 0.001, etc.
Si resolvemos las derivadas parciales obtenemos lo siguiente.
\[a=a-\alpha \frac{1}{m}\sum_{i=1}^{m}(h(x^{i})-y^{i})\]
\[b=b-\alpha \frac{1}{m}\sum_{i=1}^{m}((h(x^{i})-y^{i})x^{i})\]
En teoría, cada vez que recalculamos estas dos ecuaciones, estaremos más y más cerca del mínimo de la función de error, y por lo tanto, nuestra recta se ajustará mejor a la nube de puntos.
He hecho una implementación en Python del algoritmo tal y como os lo he contado, para que sea más claro. Sin embargo, este tipo de implementación no es muy eficiente. Normalmente se usa cálculo matricial, pero me ha parecido más pedagógico aplicar directamente las ecuaciones que acabamos de ver. Empecemos con la función para calcular el error (también llamada función de coste o loss function).

def coste(x, y, a, b):
    m = len(x)
    error = 0.0
    for i in range(m):
        hipotesis = a+b*x[i]
        error +=  (y[i] - hipotesis) ** 2
    return error / (2*m)

La función obtiene como parámetros una lista de valores x (metros cuadrados de las casas) e y (precio de las casas), así como los valores a y b para la hipótesis. A partir de ahí simplemente se calcula la ecuación de error que hemos visto más arriba y se devuelve el resultado.

def descenso_gradiente(x, y, a, b, alpha, epochs):
    m = len(x)
    hist_coste = []
    for ep in range(epochs):
        b_deriv = 0
        a_deriv = 0
        for i in range(m):
            hipotesis = a+b*x[i]
            a_deriv += hipotesis - y[i]
            b_deriv += (hipotesis - y[i]) * x[i]
            hist_coste.append(coste(x, y, a, b))
        a -= (a_deriv / m) * alpha
        b -= (b_deriv / m) * alpha
        
    return a, b, hist_coste

La función de descenso de gradiente, además de los parámetros de la función de coste, recibe alpha (learning rate) y el número de iteraciones, o lo que es lo mismo, el número de veces que vamos a recalcular y aplicar la pendiente de a y b para ir acercándonos al mínimo de la función de error. La función devuelve los valores óptimos de a y b que mejor ajustan la recta a la nube de puntos (o al menos un valor cercano) y una lista con el histórico de la función de coste, para poder estudiar cómo ha ido decreciendo el error. Vamos a probar nuestro algoritmo.

a=1
b=1
alpha = 0.0001
iters = 100000
a,b, hist_coste = descenso_gradiente(x, y, a, b, alpha, iters)
a,b 

(43.41993518666135, 1.08466250771225)

Por lo tanto, nuestra recta será \(y=43.41+1.08x\), que como vemos en la siguiente figura, ajusta bastante bien.

Si representamos gráficamente el error que hemos ido cometiendo en cada iteración podemos ver cómo va decreciendo (en este caso he omitido las primeras iteraciones porque al principio decrece muy bruscamente y la escala de la gráfica no permite observar bien la curva).

Así pues, si queremos saber cuanto deberíamos pagar por una casa de, digamos 100 metros, sólo hay que calcularlo sustituyendo b en la ecuación de la recta: \(y=43.41+1.08 \cdot 100=151.41\)
Es importante elegir bien el número de iteraciones del algoritmo. Si nos quedamos cortos nuestra recta no se ajustará bien, y si nos pasamos llega un momento en el que la mejora no es apreciable. Un método posible es detener las iteraciones cuando se compruebe que la función ya no se minimiza más allá de un umbral pequeño en cada iteración. Para que os hagáis una idea de como afecta este parámetro os dejo una animación en la que se muestra la recta calculada para diferentes valores del número de iteraciones.

Este artículo trata de ser didáctico, por lo que si quieres usar regresión lineal en la vida real, te aconsejo que mejor uses alguna buena librería como Scikit-learn que tiene un estupendo soporte para hacer este tipo de regresión.
Os dejo el enlace a mi GitHub con el notepad de Jupyter que contiene todo el código de este artículo: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/Regresion%20lineal.ipynb