El Proyecto | Sobre nosotros | Contribuir | Donaciones | Licencia

HOME

En 05-02 clasificamos: retraso sí o no. Pero el jefe de operaciones de Rutalia quiere más: "no me digas si llegará tarde, dime cuántos minutos tardará". Predecir un valor continuo es regresión, la otra mitad del aprendizaje supervisado. En esta lección la trabajamos a fondo con el dataset canónico: la regresión lineal con su solución analítica exacta (álgebra pura con numpy), y después el protagonista absoluto, el descenso de gradiente — el algoritmo iterativo de optimización que quedó pendiente en 05-02 y que es, literalmente, el motor de las redes neuronales de 05-04. Completaremos el cuadro con regularización Ridge/Lasso, regresión polinómica (donde el trade-off sesgo-varianza de 05-01 se formaliza) y las métricas propias de la regresión.

Contenido

El problema: predecir minutos_entrega
Regresión lineal simple: la recta de mínimos cuadrados
Regresión lineal múltiple: solución analítica con numpy
Descenso de gradiente: EL algoritmo
Batch, estocástico y mini-batch
Regularización: Ridge y Lasso
Regresión polinómica y el trade-off sesgo-varianza
Métricas de regresión: MSE, RMSE, MAE, R²

El problema: predecir minutos_entrega

Mismo dataset canónico de 05-01 (semilla 42), nuevo objetivo: la columna continua minutos_entrega.

datos = generar_dataset()                       # 05-01
X_num = np.column_stack([datos["distancia_km"], datos["peso_kg"],
                         datos["hora_salida"], datos["dia_semana"]])
X = np.column_stack([X_num, one_hot(datos["zona"], ZONAS)])
y = datos["minutos_entrega"]                    # ← continuo, ya no binario
X_tr, X_te, y_tr, y_te = train_test_split_manual(X, y)
X_tr_e, X_te_e = estandarizar(X_tr, X_te)

Recuerda que la verdad oculta del generador es multiplicativa (base × congestión × hora punta ×...): un modelo lineal no podrá capturarla del todo. Perfecto: veremos con precisión cuánto se le escapa y por qué.

Regresión lineal simple: la recta de mínimos cuadrados

Empecemos con una sola feature, la más informativa: distancia_km. Buscamos la recta

minutos ≈ w · distancia + b

¿Qué significa "la mejor recta"? La intuición geométrica: dibuja las 2000 entregas como puntos (distancia, minutos). Cada recta candidata deja un residuo en cada punto: la distancia vertical entre el minuto real y el predicho. Mínimos cuadrados elige la recta que minimiza la suma de los residuos al cuadrado. ¿Por qué al cuadrado y no el valor absoluto? Tres razones: penaliza más los errores grandes, hace la función derivable en todas partes (clave para lo que viene) y tiene solución cerrada:

w = cov(x, y) / var(x)        b = ȳ − w·x̄

x = datos["distancia_km"]
w = np.cov(x, y, bias=True)[0, 1] / x.var()   # pendiente
b = y.mean() - w * x.mean()                    # ordenada en el origen
print(f"minutos ≈ {w:.2f}·km + {b:.2f}")

Obtendrás algo como minutos ≈ 6.9·km + 6.5: el modelo ha redescubierto aproximadamente los "6 min/km + 5 fijos" del generador — sin que nadie se lo dijera. La pendiente no es exactamente 6.0 porque los factores multiplicativos (congestión, hora punta) inflan la relación media. Esto es aprender de los datos en su forma más pura.

Supuestos que conviene conocer (el modelo lineal los asume): relación realmente lineal, residuos independientes, de varianza constante y sin outliers dominantes. Cuando fallan, la recta sigue calculándose... pero deja de significar lo que crees.

Regresión lineal múltiple: solución analítica con numpy

Con las 13 features, buscamos el vector w (13 pesos) y el sesgo b:

ŷ = X·w + b

El truco estándar es absorber b añadiendo una columna de unos a X (así b es un peso más). Minimizar la suma de cuadrados ‖X·w − y‖² tiene solución analítica, las ecuaciones normales:

w* = (XᵀX)⁻¹ Xᵀ y

La intuición geométrica en una frase: X·w solo puede producir vectores en el subespacio generado por las columnas de X; el mejor ŷ es la proyección ortogonal de y sobre ese subespacio, y las ecuaciones normales dicen exactamente eso (el residuo debe ser perpendicular a todas las columnas: Xᵀ(y − Xw) = 0).

def regresion_lineal_analitica(X, y):
    """Resuelve las ecuaciones normales. Devuelve pesos (incluye el sesgo)."""
    Xb = np.column_stack([np.ones(len(X)), X])     # columna de unos → sesgo
    # lstsq es más estable numéricamente que invertir XᵀX explícitamente
    w, *_ = np.linalg.lstsq(Xb, y, rcond=None)
    return w

def predecir_lineal(w, X):
    return np.column_stack([np.ones(len(X)), X]) @ w

w = regresion_lineal_analitica(X_tr_e, y_tr)
pred = predecir_lineal(w, X_te_e)
print(f"RMSE test: {np.sqrt(np.mean((pred - y_te) ** 2)):.2f} minutos")

Coste: O(n·d²) para formar XᵀX más O(d³) para resolver. Con d=13 es instantáneo. Entonces, ¿para qué complicarse con métodos iterativos? Porque la solución cerrada no escala ni generaliza: con d = 10⁵ features el O(d³) explota, y —más importante— solo existe para mínimos cuadrados lineales. Para la regresión logística de 05-02 no hay fórmula cerrada. Para las redes de 05-04, tampoco. Necesitamos un algoritmo general de minimización.

Descenso de gradiente: EL algoritmo

Este apartado es el corazón de la lección — y el puente hacia todo lo que queda de módulo.

La función de coste

Definimos el error cuadrático medio (MSE) como función de los parámetros:

J(w, b) = (1/n) · Σᵢ (ŷᵢ − yᵢ)²     con  ŷᵢ = w·xᵢ + b

Cambia el punto de vista: los datos están fijos; las variables son los parámetros. J es una superficie sobre el espacio de parámetros — para mínimos cuadrados, un paraboloide convexo con un único mínimo (el que la fórmula analítica encuentra de golpe).

La idea: bajar la montaña a ciegas

Imagina estar en esa superficie con niebla: solo percibes la pendiente bajo tus pies. Estrategia: da un paso pequeño cuesta abajo y repite. El gradiente ∇J (vector de derivadas parciales) apunta en la dirección de máxima subida; caminar en −∇J es la bajada más rápida localmente.

repetir hasta converger:
    w ← w − α · ∂J/∂w
    b ← b − α · ∂J/∂b

α es el learning rate (tasa de aprendizaje): el tamaño del paso.

Las derivadas

Derivando J (regla de la cadena sobre el cuadrado; el 2 se absorbe en la convención):

∂J/∂w = (2/n) · Xᵀ(ŷ − y)          ∂J/∂b = (2/n) · Σ(ŷ − y)

Lectura intuitiva: el gradiente respecto a cada peso es el error medio ponderado por la feature correspondiente. Si el modelo se queda corto en las entregas largas, la componente de distancia_km del gradiente empuja su peso hacia arriba. El error, literalmente, dirige el aprendizaje.

Implementación con traza del descenso

def descenso_gradiente(X, y, alpha=0.1, epocas=200):
    """Minimiza el MSE. Devuelve pesos, sesgo e historial de coste."""
    n, d = X.shape
    w, b = np.zeros(d), 0.0                    # arrancamos en el origen
    historial = []
    for epoca in range(epocas):
        y_hat = X @ w + b                      # 1. predicción actual
        error = y_hat - y                      # 2. residuos
        grad_w = (2 / n) * (X.T @ error)       # 3. gradiente respecto a w
        grad_b = (2 / n) * error.sum()         #    ...y respecto a b
        w -= alpha * grad_w                    # 4. paso cuesta abajo
        b -= alpha * grad_b
        historial.append(np.mean(error ** 2))  # coste de esta época
    return w, b, historial

w_gd, b_gd, hist = descenso_gradiente(X_tr_e, y_tr, alpha=0.1, epocas=200)
for e in (0, 9, 49, 199):
    print(f"época {e:3d}: MSE = {hist[e]:8.2f}")

Verás el MSE desplomarse en las primeras épocas y estabilizarse: eso es la convergencia. Compara w_gd con la solución analítica: deben coincidir hasta el tercer decimal — dos algoritmos radicalmente distintos, el mismo mínimo, porque el paraboloide solo tiene uno.

El learning rate: la pieza delicada

α	Comportamiento	Síntoma en el historial de coste
Demasiado pequeño (0.0001)	Converge, pero tras muchísimas épocas	Baja con lentitud exasperante
Adecuado (≈0.01-0.3 aquí)	Convergencia rápida y estable	Caída suave hasta una meseta
Demasiado grande (1.5)	Cada paso salta al otro lado del valle, cada vez más lejos	El coste crece y diverge a infinito/NaN

Prueba alpha=1.5 y observa la explosión: es el error más común en la práctica. El historial de coste es tu electrocardiograma: míralo siempre. Nota también por qué estandarizamos: con features en escalas dispares, la superficie J es un valle alargado y estrecho, y el descenso zigzaguea; con features estandarizadas el valle es redondo y el mismo α sirve para todas las direcciones.

Por qué este algoritmo es EL algoritmo

El descenso de gradiente solo exige que el coste sea derivable. Cambia el MSE por el coste de la regresión logística: mismas cuatro líneas. Encadena millones de parámetros en una red neuronal: el mismo bucle, con el gradiente calculado por backpropagation (05-04). Todo el deep learning moderno es, en esencia, este for con gradientes calculados con astucia. Por eso lo hemos escrito a mano.

Batch, estocástico y mini-batch

Nuestra versión usa las n entregas en cada paso (batch). Con los millones de registros reales de Rutalia, cada paso costaría demasiado. Alternativas:

Variante	Datos por paso	Coste por paso	Trayectoria	Uso típico
Batch	Todos (n)	O(n·d)	Suave, determinista	Datasets pequeños
Estocástico (SGD)	1 ejemplo	O(d)	Muy ruidosa, oscila alrededor del mínimo	Streaming, online
Mini-batch	32-256 ejemplos	O(B·d)	Ruido moderado	El estándar (y el de 05-04)

El ruido del SGD no es solo un defecto: en superficies no convexas (las de las redes neuronales) ayuda a escapar de mínimos locales pobres. El mini-batch equilibra: gradiente razonablemente fiable, coste por paso acotado y operaciones matriciales que aprovechan la vectorización. Guarda este concepto: reaparece tal cual en el entrenamiento de la red de 05-04.

def descenso_minibatch(X, y, alpha=0.05, epocas=50, batch=64, semilla=42):
    rng = np.random.default_rng(semilla)
    n, d = X.shape
    w, b = np.zeros(d), 0.0
    for _ in range(epocas):
        orden = rng.permutation(n)                 # barajar cada época
        for i in range(0, n, batch):
            lote = orden[i:i + batch]
            err = X[lote] @ w + b - y[lote]
            w -= alpha * (2 / len(lote)) * (X[lote].T @ err)
            b -= alpha * (2 / len(lote)) * err.sum()
    return w, b

Regularización: Ridge y Lasso

Con 13 features vamos sobrados de datos. Pero imagina el caso real: cientos de features (una por cruce de zona×hora×día, meteorología, tipo de cliente...). Con muchas features y datos limitados, el modelo lineal encuentra pesos enormes que explotan correlaciones espurias del entrenamiento: sobreajuste, otra vez.

La regularización añade al coste una penalización por pesos grandes:

Método	Coste	Efecto sobre los pesos
Ridge (L2)	MSE + λ·Σwⱼ²	Los encoge todos suavemente hacia 0; nunca llegan a 0 exacto
Lasso (L1)	MSE + λ·Σ\|wⱼ\|	Lleva pesos irrelevantes a exactamente 0: selecciona features

λ controla la fuerza: λ=0 es la regresión normal; λ enorme aplasta todo hacia el modelo trivial. Se elige con validación cruzada (05-01), nunca con el test.

¿Por qué L1 produce ceros exactos y L2 no? La penalización L2 es un paraboloide (su gradiente se anula suavemente en 0: empuja proporcionalmente al peso, cada vez menos), mientras que L1 tiene "esquinas": empuja hacia 0 con fuerza constante λ, gane lo que gane el peso, y los pesos que aportan menos que λ acaban clavados en cero. Por eso Lasso sirve como selector automático de features: en Rutalia te diría, con cientos de columnas candidatas, cuáles importan de verdad.

En Ridge, incorporarlo al descenso de gradiente cuesta una línea (grad_w += 2 * lam * w), y hasta conserva solución analítica: w = (XᵀX + λI)⁻¹Xᵀy. En sklearn: Ridge(alpha=...), Lasso(alpha=...) (ojo: sklearn llama alpha a λ, no al learning rate).

Regresión polinómica y el trade-off sesgo-varianza

La verdad de Rutalia no es lineal (los factores se multiplican). Truco elegante: mantener el modelo lineal en los pesos pero enriquecer las features con términos derivados: distancia², distancia·hora_punta, etc. La maquinaria (analítica o gradiente) no cambia en absoluto: solo crece X.

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

for grado in (1, 2, 5):
    modelo = make_pipeline(PolynomialFeatures(grado), LinearRegression())
    modelo.fit(X_tr_e, y_tr)
    mse_tr = np.mean((modelo.predict(X_tr_e) - y_tr) ** 2)
    mse_te = np.mean((modelo.predict(X_te_e) - y_te) ** 2)
    print(f"grado {grado}: MSE train={mse_tr:7.2f}  test={mse_te:7.2f}")

Grado 1 se queda corto; grado 2 captura las interacciones multiplicativas del generador y mejora el test; grado 5 clava el train y empeora el test. Es la curva en U de 05-01, ahora con nombre técnico. El error de generalización se descompone en:

Sesgo: error por rigidez del modelo — la recta no puede representar la verdad multiplicativa (grado 1: sesgo alto).
Varianza: error por sensibilidad a la muestra concreta — el polinomio de grado 5 cambiaría mucho si regenerásemos los datos con otra semilla (varianza alta).
Ruido irreducible: el rng.normal(0, 4) del generador. Ningún modelo bajará de un MSE ≈ 16 (= 4²) en test; si alguno lo hace en train, está memorizando ruido.

Capacidad ↑ ⇒ sesgo ↓ y varianza ↑: el óptimo está en el codo. Todo lo que hemos visto para combatir el sobreajuste (k de k-NN, poda de árboles, bagging, λ de Ridge/Lasso, grado del polinomio) son perillas del mismo dial.

Mención breve: los árboles de regresión son los árboles de 05-02 con hojas que predicen la media de sus ejemplos y cortes que minimizan varianza en vez de Gini; el random forest de regresión promedia sus salidas. Capturan no linealidades sin ingeniería de features y son la alternativa fuerte al enfoque polinómico.

Métricas de regresión: MSE, RMSE, MAE, R²

Métrica	Fórmula	Unidades	Lectura
MSE	media de (ŷ−y)²	minutos²	La que optimizamos; penaliza mucho los errores grandes
RMSE	√MSE	minutos	La misma información, en unidades interpretables
MAE	media de \|ŷ−y\|	minutos	Error "típico"; robusta ante outliers
R²	1 − MSE/var(y)	adimensional	Fracción de la varianza explicada; 1=perfecto, 0=como predecir la media, <0=peor que la media

def metricas_regresion(y_real, y_pred):
    mse = np.mean((y_pred - y_real) ** 2)
    return {"MSE": mse, "RMSE": np.sqrt(mse),
            "MAE": np.mean(np.abs(y_pred - y_real)),
            "R2": 1 - mse / np.var(y_real)}

Para hablar con negocio, usa RMSE o MAE ("nos equivocamos unos 5 minutos por entrega"); para comparar modelos entre sí, R². Si RMSE ≫ MAE, hay outliers: unas pocas entregas con errores enormes que el cuadrado amplifica — investígalas antes de cambiar de modelo.

Errores Comunes y Consejos

Learning rate a ciegas. Si el coste diverge, α es demasiado grande; si baja a paso de tortuga, demasiado pequeño. Grafica el historial de coste en cada experimento: es tu herramienta de diagnóstico número uno.
Descenso de gradiente sin estandarizar. Con escalas dispares, el valle de J es estrecho y alargado y un α único no sirve para todas las direcciones: convergencia lentísima o divergencia. Estandariza siempre antes de gradiente.
Regularizar el sesgo b. La penalización debe aplicarse solo a los pesos w; penalizar b desplaza artificialmente las predicciones. (sklearn ya lo hace bien; en tu implementación manual, cuida ese detalle.)
Elegir λ o el grado del polinomio mirando el test. Es la fuga de 05-01 en su disfraz más frecuente. Hiperparámetros → validación cruzada; test → una sola vez, al final.
Reportar MSE a negocio. "MSE de 26 minutos cuadrados" no significa nada para operaciones. Traduce a RMSE/MAE en minutos.
Consejo: ten siempre dos referencias: el modelo trivial (predecir la media, R²=0) por abajo y el ruido irreducible por abajo del error (aquí RMSE ≈ 4): entre ambos vive todo modelo honesto.

Ejercicios

El electrocardiograma del descenso. Ejecuta descenso_gradiente sobre el dataset canónico estandarizado con alpha ∈ {0.001, 0.1, 1.1} y 100 épocas. Imprime el MSE en las épocas 0, 10, 50 y 99 para cada α. Clasifica cada comportamiento (lento / correcto / divergente) y comprueba que con α=0.1 los pesos coinciden con la solución analítica.
Ridge contra el sobreajuste polinómico. Genera features polinómicas de grado 5 y compara LinearRegression() con Ridge(alpha=10) sobre train y test (MSE). ¿Cuánto reduce Ridge el hueco train-test? Prueba también Lasso(alpha=0.5) y cuenta cuántos coeficientes deja exactamente en 0 (np.sum(modelo.coef_ == 0)).
¿Cuánto margen queda? Calcula RMSE y R² en test para: (a) el modelo trivial que predice la media de train, (b) tu regresión lineal, (c) un RandomForestRegressor(n_estimators=200). Sabiendo que el ruido del generador es σ=4, ¿qué RMSE es inalcanzable y cuánto le queda por exprimir a cada modelo?

Soluciones

Ejercicio 1:

for a in (0.001, 0.1, 1.1):
    w_, b_, h = descenso_gradiente(X_tr_e, y_tr, alpha=a, epocas=100)
    print(f"α={a}: " + "  ".join(f"e{e}={h[e]:.1f}" for e in (0, 10, 50, 99)))

w_analitica = regresion_lineal_analitica(X_tr_e, y_tr)
w_gd, b_gd, _ = descenso_gradiente(X_tr_e, y_tr, alpha=0.1, epocas=2000)
print(np.allclose(w_analitica[1:], w_gd, atol=1e-2))   # True

α=0.001: el MSE baja pero en la época 99 sigue lejos del mínimo (lento). α=0.1: caída rápida y meseta (correcto). α=1.1: el MSE crece época a época hasta overflow (divergente: cada paso salta más allá del otro lado del valle). Con suficientes épocas, α=0.1 reproduce la solución analítica: dos caminos, un solo mínimo, porque el MSE lineal es convexo.

Ejercicio 2:

from sklearn.linear_model import Ridge, Lasso

pf = PolynomialFeatures(5)
Xp_tr, Xp_te = pf.fit_transform(X_tr_e), pf.transform(X_te_e)

for nombre, m in [("lineal", LinearRegression()),
                  ("ridge", Ridge(alpha=10)), ("lasso", Lasso(alpha=0.5))]:
    m.fit(Xp_tr, y_tr)
    tr = np.mean((m.predict(Xp_tr) - y_tr) ** 2)
    te = np.mean((m.predict(Xp_te) - y_te) ** 2)
    print(f"{nombre}: train={tr:.1f} test={te:.1f}")
print("coefs a cero en lasso:", np.sum(Lasso(alpha=0.5).fit(Xp_tr, y_tr).coef_ == 0))

La lineal de grado 5 muestra un hueco train-test notable (sobreajuste); Ridge encoge los pesos y cierra buena parte del hueco cediendo un poco de train. Lasso, además, anula la gran mayoría de los cientos de coeficientes polinómicos: los términos supervivientes señalan las interacciones que importan de verdad (distancia×hora, distancia×zona congestionada).

Ejercicio 3:

from sklearn.ensemble import RandomForestRegressor

candidatos = {
    "trivial": np.full(len(y_te), y_tr.mean()),
    "lineal": predecir_lineal(regresion_lineal_analitica(X_tr_e, y_tr), X_te_e),
    "bosque": RandomForestRegressor(n_estimators=200, random_state=42)
              .fit(X_tr, y_tr).predict(X_te),
}
for nombre, pred in candidatos.items():
    m = metricas_regresion(y_te, pred)
    print(f"{nombre}: RMSE={m['RMSE']:.2f}  R2={m['R2']:.3f}")

El trivial marca el suelo (R²≈0, RMSE = desviación típica de y). Ningún modelo puede bajar de RMSE ≈ 4 (el σ del ruido del generador: es información que no existe en las features). La lineal se queda a medio camino (sesgo: no representa la verdad multiplicativa); el bosque, que sí captura no linealidades, se acerca mucho más al suelo de 4. La distancia de cada RMSE a 4 es el margen de mejora real de cada modelo.

Conclusión

Esta lección deja tres capas bien asentadas. La regresión lineal con su doble solución: la analítica (proyección ortogonal, ecuaciones normales, lstsq) y la iterativa. El descenso de gradiente desarrollado pieza a pieza —coste MSE, derivadas, learning rate, historial de convergencia, variantes batch/estocástico/mini-batch— y verificado contra la solución exacta. Y el control del sobreajuste con vocabulario formal: sesgo contra varianza, con la regularización Ridge/Lasso y el grado polinómico como perillas, y el ruido irreducible como suelo que ningún modelo honesto perfora. Pero fíjate en la limitación que arrastramos: para capturar la verdad no lineal de Rutalia hemos tenido que fabricar las features nosotros (polinomios, interacciones) — otra vez conocimiento aportado a mano. ¿Y si el modelo pudiera aprender también las transformaciones? Eso es exactamente una red neuronal: capas de regresiones (lineales + activaciones) apiladas, entrenadas de extremo a extremo con el mismo descenso de gradiente que acabas de implementar. El gradiente, eso sí, habrá que propagarlo hacia atrás a través de las capas: backpropagation. Nos vemos en 05-04.

Algoritmos de Regresión