Proyecciones: lo que hace la regresión lineal por dentro

La sombra de un vector sobre otro 🌸

Imagínate una linterna justo encima de un palo apoyado en la pared. La sombra que cae sobre la pared es más corta que el palo, y va en la dirección de la pared. Eso es una proyección.

En vectores es lo mismo: dado a, ¿cuánto de a va en la dirección de b?

Volvemos a abrir el archivo

Cada capítulo arranca de cero, así que estas líneas se repiten. Son las mismas del capítulo 1 y las vas a ver en todos.

import pandas as pd
import numpy as np

df = pd.read_csv('ventas-miss-yera.csv')
for c in ('unidades', 'monto', 'satisfaccion'):
    df[c] = pd.to_numeric(df[c], errors='coerce')
X = df[['unidades', 'monto', 'satisfaccion']].dropna().to_numpy()
print('X.shape:', X.shape)

X.shape: (2238, 3)

Proyectar, que son dos productos punto

a, b = X[0], X[1]
coef = float(a @ b) / float(b @ b)
proy = coef * b

print('a          :', np.round(a, 2))
print('b          :', np.round(b, 2))
print('coeficiente:', round(coef, 6))
print('proyeccion :', np.round(proy, 2))

a          : [ 10.   480.37   4.  ]
b          : [ 10.  524.9   3. ]
coeficiente: 0.915209
proyeccion : [  9.15 480.39   2.75]

Ese 0.915 quiere decir "a es como el 91.5% de b, en esa dirección". Y la proyección es b estirado por ese número, o sea la operación de escalar del capítulo 3.

Fíjate en que la proyección va en la dirección de b. No es un pedazo de a: es b encogido para que se parezca lo más posible a a.

Lo que sobra queda perpendicular

resto = a - proy
print('resto     :', np.round(resto, 2))
print('resto . b :', round(float(resto @ b), 8))

resto     : [ 0.85 -0.02  1.25]
resto . b : -0.0

Producto punto cero quiere decir perpendicular, y no es casualidad: es la definición. Dos vectores son ortogonales cuando su producto punto da cero, punto.

Y esto pasa siempre, no solo con estos números. La proyección es, por construcción, la parte que sí se puede explicar con b. Lo que queda es justamente lo que b no puede explicar, y por eso no le queda nada de b dentro 🐣

Y aquí está la regresión

Hasta ahora proyectamos sobre un vector. Ahora proyecto sobre todas las columnas a la vez: quiero predecir el monto con las unidades y la satisfacción.

y = X[:, 1]
A = np.column_stack([np.ones(X.shape[0]), X[:, 0], X[:, 2]])
w = np.linalg.solve(A.T @ A, A.T @ y)

print('intercepto  :', round(float(w[0]), 4))
print('unidades    :', round(float(w[1]), 4))
print('satisfaccion:', round(float(w[2]), 4))

intercepto  : 833.6427
unidades    : 0.463
satisfaccion: -11.819

Esa columna de unos del principio es el intercepto, la b de la fórmula del capítulo 4. Meterla como una columna más es el truco para que la fórmula la trate igual que a las otras.

Y ahora la comprobación que quiero que veas:

yhat = A @ w
res = y - yhat
print('yhat[:3] :', np.round(yhat[:3], 2))
print('residuos contra cada columna:', np.round(A.T @ res, 6))

yhat[:3] : [791.   802.82 827.84]
residuos contra cada columna: [ 0. -0.  0.]

Ceros. Los residuos son ortogonales a las tres columnas, y eso es exactamente lo mismo que pasó arriba con dos vectores, solo que en tres direcciones a la vez.

Eso es lo que significa "mínimos cuadrados": de todos los puntos que se pueden construir con tus columnas, la predicción es el más cercano a la realidad. Y el más cercano es siempre la sombra 🌸

Ahora la parte incómoda

print('RMSE del modelo         :',
      round(float(np.sqrt((res ** 2).mean())), 2))
print('RMSE de predecir la media:',
      round(float(np.sqrt(((y - y.mean()) ** 2).mean())), 2))

RMSE del modelo         : 753.47
RMSE de predecir la media: 753.67

753.47 contra 753.67. El modelo mejoró veinte céntimos sobre no hacer nada.

Y quiero que te quedes con esto porque es lo que un curso normal te esconde: la proyección funcionó perfecto y el modelo no sirve. Las dos cosas a la vez.

El álgebra hizo su trabajo impecablemente: encontró la mejor combinación posible de esas columnas. Lo que pasa es que en las unidades y la satisfacción de este archivo no hay información sobre el monto, así que la mejor combinación posible sigue siendo mala.

Los residuos ortogonales son la firma de que ajustaste bien. No son una nota del modelo. Confundir las dos cosas es de las equivocaciones más caras que hay en este oficio 🙃

Lo que te llevas

• Proyectar es quedarse con el trozo que va en una dirección, y son dos productos punto.
• Producto punto cero quiere decir perpendicular.
• Lo que sobra después de proyectar siempre queda perpendicular.
• Una regresión por mínimos cuadrados es una proyección sobre tus columnas.
• La columna de unos es el intercepto.
• Residuos ortogonales dicen que ajustaste bien, no que el modelo sirva.

Comprueba que se entendió

Ejercicios

a, b = X[0], X[1]
proy = (float(a @ b) / float(b @ b)) * b
resto = a - proy
print('|a|^2            :', round(float(a @ a), 4))
print('|proy|^2 + |resto|^2:',
      round(float(proy @ proy + resto @ resto), 4))

|a|^2            : 230871.3369
|proy|^2 + |resto|^2: 230871.3369

y = X[:, 1]
u = X[:, 0]
coef = float(u @ y) / float(u @ u)
print('coeficiente:', round(coef, 4))
print('o sea que cada unidad vale', round(coef, 2), 'soles segun este ajuste')
res = y - coef * u
print('residuos contra u:', round(float(u @ res), 6))

coeficiente: 55.3417
o sea que cada unidad vale 55.34 soles segun este ajuste
residuos contra u: 0.0

rng = np.random.default_rng(0)
basura = rng.normal(size=X.shape[0])
y = X[:, 1]
A = np.column_stack([np.ones(X.shape[0]), basura])
w = np.linalg.solve(A.T @ A, A.T @ y)
res = y - A @ w
print('residuos contra cada columna:', np.round(A.T @ res, 6))
print('RMSE :', round(float(np.sqrt((res ** 2).mean())), 2))

residuos contra cada columna: [0. 0.]
RMSE : 753.67