Fórmulas y implementación desde cero de regresión lineal simple y múltiple

Las fórmulas en sí son muy cortas. Como primera lección de ESL, es muy básico.

Hay muchas bibliotecas que pueden hacer esto, pero recientemente necesité calcular regresión con GPU en mi desarrollo (lo cual puede acelerar muchas veces), así que aproveché para escribir un blog y probar si la visualización de fórmulas del sitio web funciona correctamente.

¿Qué es la regresión lineal?

En pocas palabras, es usar una ecuación lineal para representar la tendencia de los datos.

Univariada: $y=mx+b$ Multivariada: $y=m_1x_1+m_2x_2+…+b$

|1.5x Puedes imaginarlo como si cada punto fuera una estrella de la misma masa, y la línea de regresión es una vara larga colocada en equilibrio final en este entorno gravitacional.

Fórmula de regresión lineal simple univariada (Linear regression)

La regresión Simple OLS se presenta por separado porque es más intuitiva y fácil de entender:

m=\frac {\sigma_{x, y}} {\sigma^2_x}=\frac {E[(X-\mu)(Y-\nu)]} { E[(X-\mu)^2] }

Es simplemente la covarianza dividida por la varianza de x. No hablaremos de cómo derivar esta fórmula, pero ya debes sentir el encanto de esta fórmula.

Luego, b se obtiene sustituyendo m:

b=\bar{y} - m\bar{x}

Implementación del código con pyTorch

Usamos pyTorch para la aceleración por GPU. Supongamos que tenemos un conjunto de datos y con los precios de NVDA de febrero de 2018:

import torch
y = torch.tensor([
    225.58, 228.8 , 217.52, 230.93, 228.03, 232.63, 241.42, 246.5 ,
    243.84, 249.08, 241.51, 242.15, 245.93, 246.58, 246.06, 242.  ,
    232.21, 236.54, 235.65, 242.16
])

Luego generamos x:

x = torch.arange(len(y)).float()

tensor([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19.])

Código de regresión (aquí solo para legibilidad, hay cálculos repetidos):

demean_x = x - x.mean()
demean_y = y - y.mean()
n_1 = x.shape[0] - 1
m = torch.sum(demean_x * demean_y / n_1) / torch.sum(demean_x ** 2 / n_1)
b = y.mean() - m * x.mean()
print(m, b)

(tensor(0.7734), tensor(230.4090))

Grafiquemos para ver el resultado:

from matplotlib import pyplot as plt
plt.plot(x, y)
plt.plot(x, m * x + b)
plt.show()

|1.5x Se ve bien 🥂.

Fórmula de regresión lineal múltiple (Multiple Linear regression)

Es esencialmente lo mismo que la regresión simple. Aquí usamos notación matricial, que es aún más concisa:

b,m_1,...,m_n=(X^TX)^{-1}X^Ty

Una sola fórmula lo resuelve todo, y también es aplicable a la regresión lineal simple.

Fuera de tema, $X^TX$ aparece en todas partes, e incluso hay un fondo de capital privado llamado XTX.

Implementación del código con pyTorch

La fórmula contiene mucha información oculta, y las explicaciones textuales son demasiado tediosas. Implementémoslo con código.

Siempre he pensado que si las fórmulas en libros de texto y artículos de investigación tuvieran el código correspondiente y los resultados de cada paso, habría muchas menos dudas al leer, ya que para que el código sea ejecutable se necesita información completa.

Aquí usamos regresión lineal bivariada para demostrar. Lo más común es usar una ecuación cuadrática univariada ($y=ax^2+bx+c$), es decir, ajustar una curva. Si es una ecuación univariada, ¿por qué es regresión bivariada? Porque los dos términos en la regresión bivariada se refieren a 2 variables independientes: $x, x^2$

Todavía usando los datos x, y anteriores, primero generemos la matriz X, que contiene $x, x^2$, pero ten en cuenta que aquí la primera columna se añade una constante 1, para que la matriz sea definida positiva y para calcular b (en la imagen ${\hat\beta}_0$):

X = torch.stack([torch.ones(x.shape), x, x ** 2]).T

tensor([[ 1., 0., 0.], [ 1., 1., 1.], [ 1., 2., 4.], [ 1., 3., 9.], [ 1., 4., 16.], [ 1., 5., 25.], … [ 1., 15., 225.], [ 1., 16., 256.], [ 1., 17., 289.], [ 1., 18., 324.], [ 1., 19., 361.]])

Luego simplemente calculamos:

b, m1, m2 = (X.T @ X).inverse() @ X.T @ y
print(b, m1, m2)

(tensor(220.5210), tensor(4.0694), tensor(-0.1735))

Sí, este es el resultado. Grafiquemos para ver:

plt.plot(x, y)
plt.plot(x, m1 * x + m2 * x**2 + b)
plt.show()

|1.5x

Perfecto 👏, ¿verdad que es simple? Pero en realidad hay conceptos más importantes detrás, como el caso de múltiples soluciones y el procesamiento ortogonal.

Colinealidad y problemas de ortogonalidad

En el uso práctico, no es común que los problemas se resuelvan simplemente con regresión polinómica. Generalmente se hace regresión sobre valores observados, y estos valores observados básicamente no son completamente independientes, es decir, la covarianza entre las x no es 0. Mantener las x independientes permite que el coeficiente sea responsable solo de esa x y no se vea afectado por otras x, lo que requiere un ajuste ortogonal.

No entraré en detalles. Si quieres resolver el problema rápidamente, aquí hay una solución aproximada: usa descomposición QR. Calculando Q, R = torch.qr(X) y luego usando la fórmula $b, m_1,…,m_n=R^{-1}Q^Ty$ para resolverlo.

Si estás interesado, puedes leer ESL 3.2.3 Multiple Regression from Simple Univariate Regression, que explica la regresión múltiple a partir de la regresión univariada simple.

Autor: Zhang Jianhao Heerozh (heerozh.com)