Método de mínimos cuadrados

El método de mínimos cuadrados busca los coeficientes de una función que mejor ajusta pares de observaciones (x_i, y_i), minimizando la suma de los residuos al cuadrado. Nos centramos en el caso más simple, una recta y = wx + b, y revisamos la intuición y una implementación práctica.

Las fórmulas se muestran con KaTeX. \(\hat y\) denota la predicción del modelo y \(\epsilon\) denota el ruido.

Objetivo #

Aprender la recta \(\hat y = wx + b\) que mejor se ajusta a los datos.
“Mejor” significa minimizar la suma de errores cuadráticos (SSE): \(\displaystyle L(w,b) = \sum_{i=1}^n (y_i - (w x_i + b))^2\)

Crear un conjunto de datos sencillo #

Generamos una recta con ruido y fijamos la semilla para reproducibilidad.

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # opcional para etiquetas en japonés

rng = np.random.RandomState(42)
n_samples = 200

# Recta verdadera (pendiente 0.8, intercepto 0.5) con ruido
X = np.linspace(-10, 10, n_samples)
epsilon = rng.normal(loc=0.0, scale=1.0, size=n_samples)
y = 0.8 * X + 0.5 + epsilon

# Convertir a 2D para scikit-learn: (n_samples, 1)
X_2d = X.reshape(-1, 1)

plt.figure(figsize=(10, 5))
plt.scatter(X, y, marker="x", label="observaciones", c="orange")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend()
plt.show()

png

En scikit-learn, las características son siempre un arreglo 2D: filas = muestras, columnas = variables. Use X.reshape(-1, 1) para una sola variable.

Inspeccionar el ruido #

Veamos la distribución de epsilon.

plt.figure(figsize=(10, 5))
plt.hist(epsilon, bins=30)
plt.xlabel("$\\epsilon$")
plt.ylabel("frecuencia")
plt.show()

png

Regresión lineal (mínimos cuadrados) con scikit-learn #

Usamos sklearn.linear_model.LinearRegression.

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

model = LinearRegression()  # fit_intercept=True por defecto
model.fit(X_2d, y)

print("pendiente w:", model.coef_[0])
print("intercepto b:", model.intercept_)

y_pred = model.predict(X_2d)

# Métricas
mse = mean_squared_error(y, y_pred)
r2 = r2_score(y, y_pred)
print("MSE:", mse)
print("R^2:", r2)

# Gráfico
plt.figure(figsize=(10, 5))
plt.scatter(X, y, marker="x", label="observaciones", c="orange")
plt.plot(X, y_pred, label="recta ajustada", c="C0")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend()
plt.show()

png

No es obligatorio escalar para mínimos cuadrados, pero ayuda con problemas multivariados y regularización.

Solución en forma cerrada (referencia) #

Para \(\hat y = wx + b\):

\(\displaystyle w = \frac{\operatorname{Cov}(x,y)}{\operatorname{Var}(x)}\)
\(\displaystyle b = \bar y - w,\bar x\)

Verificar con NumPy:

x_mean, y_mean = X.mean(), y.mean()
w_hat = ((X - x_mean) * (y - y_mean)).sum() / ((X - x_mean) ** 2).sum()
b_hat = y_mean - w_hat * x_mean
print(w_hat, b_hat)

Errores comunes #

Formas de arreglos: X debe ser (n_samples, n_features). Con una variable, use reshape(-1, 1).
Forma del objetivo: y puede ser (n_samples,). (n,1) también funciona; cuide el broadcasting.
Intercepto: fit_intercept=True por defecto. Si centró X e y, False está bien.
Reproducibilidad: fije la semilla con np.random.RandomState o np.random.default_rng.

Más allá (multivariado) #

Con varias características, mantenga X como (n_samples, n_features). Un pipeline combina preprocesamiento y estimador.

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X_multi = rng.normal(size=(n_samples, 2))
y_multi = 1.0 * X_multi[:, 0] - 2.0 * X_multi[:, 1] + 0.3 + rng.normal(size=n_samples)

pipe = make_pipeline(StandardScaler(), LinearRegression())
pipe.fit(X_multi, y_multi)

El código es para aprendizaje; las figuras están prerenderizadas para el sitio.