Machine learning desde cero (II)
Regresión lineal multidimensional

Repasito unidimensional

En el post anterior vimos cómo la predicción basada en una familia de datos de salida y una sola familia de datos de entrada se podía expresar matemáticamente como el ajuste de una recta a una nube de puntos en el plano.

La función objetivo de nuestro modelo consistía en encontrar los parámetros que minimizaban la recta de mejor ajuste a la nube de puntos. Es decir, minimizar la expresión:

A esta expresión la llamábamos función objetivo. A modo de ejercicio busqué el mínimo de la función objetivo, despejé para $\alpha$ y $\beta$ y llegué a la siguiente expresión matricial:

Esta fue, precisamente, la expresión que utilicé en el script de python para realizar el cálculo de $\alpha$ y $\beta$, aprovechando que la librería numpy nos permite realizar cálculos matriciales de forma muy sencilla:

...
n = len(x_data)
sum_x = np.sum(x_data)
sum_x_squared = np.sum(np.power(x_data, 2))
sum_y = np.sum(y_data)
sum_xy = np.sum(x_data.dot(y_data))

a = np.array([[n, sum_x], [sum_x, sum_x_squared]]) # matriz que invertimos
b = np.transpose([sum_y, sum_xy]) # vector

self.w = np.linalg.inv(a).dot(b) # w es nuestro vector con los valores de alpha y beta
...

Dos mejor que uno

Supongamos ahora que queremos realizar predicciones en base a dos inputs en lugar de uno solo. Siguiendo el ejemplo de las viviendas del post anterior, supongamos que queremos predecir el precio de las viviendas no sólo en función de los metros cuadrados habitables sino también de su antigüedad en años.

Dejo por aquí una nueva tabla de ejemplo con nuestros nuevos datos:

Si quisiéramos plotear esta nueva tabla, es fácil ver que, si mantenemos nuestra variable respuesta precio en la vertical, necesitaremos un nuevo eje cartesiano horizontal para la nueva variable antigüedad. Si dos ejes cartesianos constituyen un plano, dos ejes horizontales y uno vertical constituyen el espacio tridimensional de toda la vida.

Nuestra nube de puntos ahora se encontrará, en lugar de en el plano, en el espacio. ¡Veámoslo!

Nuestra vieja función objetivo contenía dos incógnitas, $\alpha$ y $\beta$, que definían la recta de mejor ajuste en nuestro modelo de predicción, a saber: el ajuste de una recta a una nube de puntos en el plano.

Ahora la mera representación de nuestros datos exige una dimensión adicional y, por tanto, nuestro modelo tendrá que ser capaz de adaptarse y capturar la tendencia de nuestra nueva nube de puntos, que ahora se encuentra en el espacio.

Siguiendo la misma filosofía de “elemento geométrico que mejor se ajusta a la nube de puntos”, ¿qué tipo de elemento geométrico lineal crees que necesitaremos capaz de ajustarse a una nube de puntos en el espacio?

¡El plano! Con el plano podemos mantener la misma filosofía de “mínimas distancias verticales” plano-punto, sin embargo va a ser necesario replantear nuestra vieja función objetivo, puesto que cuando la construimos no teníamos en cuenta la nueva dimensión que hemos introducido. ¡Presta atención a la pauta que sigue la nueva función objetivo con respecto a la vieja, porque más adelante nos ayudará a construir la función objetivo para cualquier número de dimensiones adicionales!

• Sean dos series de inputs de entrada (dos variables predictoras) ${x_1}_i$ y ${x_2}_i$
• Sea una serie de outputs de salida (variable respuesta) $y_i$
• Y sea un plano arbitrario definido como $y=\alpha+\beta_1 x_1+\beta_2 x_2$

Para cada punto del espacio definido como $({x_1}_i, {x_2}_i, y_i)$, la distancia vertical de dicho punto a nuestro plano arbitrario será:

Ahora es un momento ideal para tratar de re-expresar nuestra función distancia de forma que facilite la introducción de nuevas variables en el futuro. En el fondo se trata de un sencillo ejercicio de álgebra e imaginación:

Para facilitar las cosas renombremos $\alpha$ como $\beta_0$ y supongamos que está multiplicando al número $1$. Nuestra función quedaría del siguiente modo:

El segundo término de la resta no es más que un sumatorio, que puede re-expresarse como el producto punto (o producto escalar) entre dos vectores, en este caso:

En la literatura matemática tradicional se puede encontrar al vector de las betas representado simplemente con una $\beta$ y con el nombre de vector de estimadores, mientras que en machine learning puede encontrarse representado por la letra $w$ y se le suele llamar vector de pesos o weights. El vector de equises (con el 1 como primer elemento), por su lado, suele representarse, sencillamente, con una $x$.

Por tanto, nuestra expresión con la notación habitual en machine learning podría re-escribirse como:

Nótese que si agregáramos ahora nuevas variables predictoras, lo único que haríamos sería incrementar la dimensión de los vectores $w$ y $x$, es decir, si en el ejemplo con dos variables predictoras sabemos que se trata de vectores $3\times1$, agregando una tercera variable se convertirían en vectores $4\times1$ y así seguiríamos hasta cualquiera número de variables predictoras. En general, para $k$ variables predictoras obtendríamos vectores de dimensión $(k+1)\times1$, sin embargo nuestra expresión para la distancia sería la misma. ¡Ya estamos un pasito más cerca de la generalización!

Vayamos ahora un poquito más allá y re-expresemos también matricialmente el conjunto de distancias, es decir, en lugar conformarnos con la expresión para cada $\epsilon_i$, pongamos todas las $epsilon$ en un solo vector y re-expresemos la parte derecha de la igualdad para satisfacer la equivalencia.

Teniendo en cuenta la expresión (1), si pasáramos todo en forma de vector, a lo loco, nos quedaría algo como esto para el caso de dos variables predictoras y n observaciones (n puntos en la nube):

Todo apunta a que la magia del álgebra lineal nos permite re-expresar el vector de más a la derecha de forma mucho más elegante, ¿qué tal un producto punto entre una matriz y el vector de betas?

En la literatura matemática encontraremos esta expresión matricial de las distancias en su forma simbólica abreviada: $\epsilon$ representando al vector de distancias, $y$ representando al vector de valores de la variable respuesta, $X$ representando a la matriz con la columna de unos y las observaciones de las variables predictoras y, finalmente, $b$ para el vector de estimadores. En machine learning, por otro lado, no es raro encontrar expresado el vector de valores de la variable respuesta como $T$ de target o valores objetivo y, como ya vimos antes, el vector de estimadores o parámetros incógnita como $w$ de vector de pesos o weights, en definitiva:

Esta expresión es elegante y potente al mismo tiempo. Primeramente porque se puede emplear para cualquier número de variables predictoras, siempre y cuando tengamos en cuenta que las dimensiones de $X$ y $w$ cambiarán. Por ejemplo: en el caso de dos variables predictoras $X$ tiene dimensiones $n\times3$ y $w$ dimensiones $3\times1$. En el caso de tres variables predictoras $X$ será una matriz $n\times4$ y $w$ un vector $4\times1$. Cae de cajón, por tanto, que en el caso general, para $n$ observaciones y $k$ variables predictoras, $X$ será una matriz $n\times(k+1)$ y $w$ un vector $(k+1)\times1$.

En segundo lugar, la gran ventaja de esta expresión matricial es que nos va a permitir expresar y resolver nuestra función objetivo sin despeinarnos y, además, ¡para cualquier número de variables predictoras!

Recordemos que nuestra función objetivo no era la expresión de las distancias en sí, sino el mínimo del sumatorio de las distancias al cuadrado, es decir:

Si abrazas la magia del álgebra lineal, este sumatorio también se puede expresar como el producto punto entre dos vectores de epsilons:

Empleando nuestra maravillosa expresión (2) al producto punto entre los vectores epsilons y operando obtendremos nuestra nueva función objetivo en forma de operaciones matriciales. Nótese que estamos multiplicando vectores y matrices, no escalares, por lo que la transposición, la propiedad distributiva y las simplificaciones se han aplicar con especial cariño (siguiendo las propiedades de vectores y matrices):

Nuestra función objetivo para el caso general será, por tanto:

Ahora hemos de derivar respecto al vector de parámetros incógnita o vector de pesos $w$. Esta derivada no es ni mucho menos trivial y consiste, de forma simplificada, en ver que no existe vector de incógnitas en el primer término (su derivada será nula), que el segundo término equivale a la derivada de una expresión lineal y que la derivada del último término equivale a la derivada de una expresión cuadrática (el vector de incógnitas aparece dos veces). Se puede consultar la demostración de estas derivadas en cualquier libro medianamente exhaustivo sobre el modelo de mínimos cuadrados (pregunta en el área de comentarios si tienes curiosidad). El resultado sería el siguiente:

Igualando a cero y despejando para el vector de incógnitas:

Nótese que esta expresión (3) nos permite calcular los parámetros $\beta$ para cualquier número de variables predictoras, incluyendo nuestro ejemplo al principio de este post (dos variables predictoras) y, también, el ejemplo del primer post sobre regresión lineal (una variable predictora). El resultado $w$ será, en cada uno de los casos, un vector con dimensiones $k+1$, donde $k$ será el número de variable predictoras con que alimentemos nuestro modelo de ajuste lineal.

Escribamos, por fin, nuestro script en python para resolver el caso general de regresión lineal y resolvamos el ejemplo que puse al principio del post, es decir, calculemos el plano que mejor se ajusta a la nube de puntos en el espacio:

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
import numpy as np


class LinearRegression(object):
    def __init__(self, number_of_predictive_variables):
        self.k = number_of_predictive_variables
        self.w = np.zeros(number_of_predictive_variables + 1)

    def train(self, y_data, *x_data):
        x_matrix = np.ones((len(y_data), self.k + 1))
        for i, x_var in enumerate(x_data, start=1):
            x_matrix[:, i] = x_var

        self.w = np.linalg.inv(x_matrix.T.dot(x_matrix)).dot(x_matrix.T).dot(y_data)

    def predict(self, *x_values):
        x_values_vector = np.ones(self.k + 1)
        for i, x_value in enumerate(x_values, start=1):
            x_values_vector[i] = x_value

        return np.array(x_values_vector).dot(self.w)

    def get_weights(self):
        return self.w


k = 2
X1 = np.array([92.7, 139.4, 109.2, 114.5, 104.1, 91.8, 115.2, 139.4, 317.7, 278.7, 113.8, 144.2, 90.6, 104.1, 94.8,
               139.4, 154.6, 138.2, 127.8, 139.4, 116.7, 157.0, 169.1, 153.5, 165.1, 139.7, 170.1, 111.5])
X2 = np.array([42, 62, 40, 54, 42, 56, 51, 32, 42, 14, 32, 30, 30, 32, 46, 32, 50, 22, 17, 23, 40, 22, 50, 44, 48, 3,
               31, 30])
Y = np.array([259000, 295000, 279000, 259000, 299000, 299000, 309000, 289000, 849000, 829000, 359000, 315000, 310000,
              309000, 300000, 289000, 369000, 419000, 405000, 439000, 375000, 379000, 445000, 379000, 389000, 369000,
              458000, 410000])

lr = LinearRegression(k)
lr.train(Y, X1, X2)

print "weights (vector of betas) =", lr.get_weights()
print "price prediction for 50 m2 and 0 years =", lr.predict(50, 0)
print "price prediction for 50 m2 and 30 years =", lr.predict(50, 30)
print "price prediction for 50 m2 and 100 years =", lr.predict(50, 100)
# weights (vector of betas) = [ 83455.24240792   2494.67526436  -1438.1104518 ]
# price prediction for 50 m2 and 0 years = 208189.005626
# price prediction for 50 m2 and 30 years = 165045.692072
# price prediction for 50 m2 and 100 years = 64377.9604464

Nuestro vector de betas o pesos define ahora un plano en el espacio. Si lo ploteamos obtenemos nuestro plano de predicción del precio de una vivienda en función de sus m2 habitables y su antigüedad:

Curiosamente, a partir de dos variables predictoras (tres, veinte, mil, etc.), ya no es posible representar gráficamente el elemento geométrico que mejor se ajusta a la nube de puntos (hiperplano en un espacio de más de tres dimensiones), de hecho, ya no nos será posible ni siquiera representar la nube de puntos. En cualquier caso, gracias a nuestro modelo generalizado para la regresión lineal multidimensional, podremos resolver numéricamente nuestro vector de pesos y predecir la variable respuesta con cualquier número de variables predictoras.

Llegados a este punto ya tenemos una idea aproximada de las bases sobre las que se asientan algunos algoritmos de clasificación que veremos más adelante. Además, podemos imaginar que los ajustes no lineales seguirán lógicas parecidas, pero utilizando figuras geométricas “flexibles”, es decir, curvas y superficies curvadas cuyo objetivo será ajustarse aún mejor a las nubes de puntos en cualesquiera dimensiones, es decir, hacer predicciones aún más fiables.

Por supuesto, estos métodos no lineales plantearán nuevas cuestiones y dificultades: ¿cómo calcularemos la función curva que mejor se ajusta a la nube de puntos? ¿si es posible ajustar perfectamente a la nube de puntos, hasta qué punto es más fiable dejar ciertas distancias figura-punto y alisar la figura geométrica curvada? ¿si hemos de alisarla, hasta qué punto hemos de alisarla? Confío en que todas estas cuestiones las iremos viendo más adelante.

Por otro lado, el análisis de la bondad de ajuste es un tema muy importante de cara a evitar algunos problemas típicos de la regresión lineal, por ejemplo:

• La posibilidad de que nuestra nube de puntos sea un cúmulo de puntos sin sentido, es decir, con muy poca correlación entre las variables predictivas y la variable respuesta. Nuestro modelo arrojará resultados en cualquier caso, pero ¿cómo de fiables serán las predicciones? ¿hay algún modo de preverlo? Para tener una idea de la fiabilidad del modelo, en el caso lineal, se calcula el llamado Coeficiente de Determinación o $R^2$.
• La posibilidad de que dos o más de las variables predictivas se “expliquen” unas a otras, por ejemplo: en el caso de las viviendas la variable metros cuadrados habitables se correlacionará estrechamente con la variable metros cuadrados construidos. “Alimentar” a nuestro modelo con variables tan estrechamente correlacionadas como estas introduce importantes distorsiones en el análisis del aporte de cada una de las variables a la bondad de ajuste del modelo. A este problema o fenómeno lo conocen nuestros amigos, los estadísticos, como multicolinealidad.

Basta decir que es posible tratar de prevenir esta clase de problemas a través de múltiples tests matemáticos que, a fin de cuentas, permiten una elección inteligente de las variables predictivas y nos ayudan, en la medida de lo posible, a incrementar la bondad de ajuste de nuestro modelo lineal. Al respecto, aconsejo echar un vistazo a cualquier libro de estadística o a cualquier curso online sobre regresión lineal aplicada para estar al tanto de estas cuestiones importantes. A mí me parece bastante práctico y al grano el curso online STAT501 del Colegio de Ciencias de la Universidad Estatal de Pensilvania.

Si alguien quiere indagar en la parte estadística y rigurosa de la regresión lineal, es bueno saber que existe software estadístico muy completo que asiste con la regresión y los tests necesarios para analizar la bondad de ajuste. A mí, personalmente, me gusta el software comercial Minitab,por sencillo e intuitivo, aunque no me cabe ninguna duda de que también se puede hacer magia estadística con el más espartano y open source R + RStudio.

En cualquier caso, nosotros vamos a seguir dándole al machine learning.

Próximamente

En el próximo post sobre machine learning compartiré mis apuntes sobre algoritmos básicos de clasificación: la clasificación lineal y la regresión logística.