20  Métodos de suavización: splines

El tráfico no cabe en una recta

Considera una simulación del flujo de una avenida a lo largo de 24 horas. Un equipo de movilidad registra el flujo vehicular (vehículos por minuto) y quiere un modelo sencillo: si me das la hora, te digo cuántos vehículos esperar. En un proyecto real, su utilidad tendría que evaluarse según la decisión de operación y el costo de errar en cada hora.

Cualquiera que haya cruzado esa avenida sabe lo que va a pasar. A las 3 de la mañana está casi vacía; a las 8 llega el pico de entrada al trabajo; al mediodía baja; a las 6 de la tarde vuelve a saturarse. El efecto de la hora no es el mismo a lo largo del día: en algunas regiones de \(X\) subir una hora casi no cambia nada, y en otras media hora marca la diferencia entre fluidez y congestión.

Una recta obliga a que el efecto de \(X\) sea constante en todo el dominio, y en muchos fenómenos eso no ocurre: el consumo eléctrico según la temperatura (la V del capítulo de regresión múltiple), el precio de una vivienda según el área cuando el mercado cambia de segmento o la concentración de un medicamento según el tiempo. La relación se dobla, se aplana o tiene más de un régimen.

En los capítulos anteriores aprendiste a evaluar si un modelo generaliza, con MAE, RMSE, \(R^2\), validación cruzada y bootstrap. Esa pregunta sigue vigente, y en este capítulo estudias una forma funcional que puede representar cambios locales como los del tráfico. Todo el capítulo gira alrededor de una pregunta de diseño:

¿Cómo permites que un modelo cambie su comportamiento en distintas regiones de \(X\) sin producir una curva artificial o discontinua?

El recorrido empieza con la recta fallando a la vista, sigue con el polinomio global, la idea de partir el dominio, los nodos (knots) y el salto en los empalmes, y termina en el spline.

• Tiempo estimado de lectura: 85 minutos (unas 3 horas si ejecutas el código y respondes las preguntas)
• Nivel: Avanzado
• Conocimientos previos recomendados: la validación cruzada k-fold, la evaluación de modelos y la regresión múltiple del nivel Avanzado, y la regresión lineal del nivel Básico.
• Tecnologías mencionadas: Python 3.13, NumPy, pandas, Matplotlib y scikit-learn (PolynomialFeatures, SplineTransformer, Pipeline).

Los datos son simulados con las fórmulas y semillas declaradas en el código: un día de tráfico con dos picos y ruido normal de desviación 6.5 vehículos por minuto, en una muestra ilustrativa de 180 mediciones y una segunda muestra independiente para seleccionar y evaluar; la función de Runge; y un día con horas muestreadas de forma desigual. No proceden de clientes. Todo corre en CPU, y las versiones se imprimen en la primera celda de código.

Objetivos

Al terminar este capítulo vas a ser capaz de:

• reconocer cuándo una relación lineal es insuficiente, en el gráfico, en los residuos y en las métricas;
• explicar qué es una regresión polinómica, por qué sigue siendo lineal en los parámetros y qué problemas tiene un polinomio global de grado alto en los extremos;
• explicar qué es ajustar por tramos, qué es un knot y en qué se diferencia un polinomio por tramos de un spline;
• describir las condiciones de continuidad de un spline cúbico y por qué importan geométricamente;
• usar SplineTransformer dentro de un Pipeline, donde el spline expande \(X\) y LinearRegression estima los pesos;
• controlar la flexibilidad con el número de knots, seleccionarla con validación cruzada y comparar recta, polinomio y spline con métricas de generalización.

Ejecuta las celdas en orden: las secciones comparten variables y funciones de dibujo.

20.1 Preparación del entorno

La celda carga NumPy, pandas, Matplotlib y scikit-learn, fija la semilla y define la paleta del libro, para que los gráficos se lean igual en todos los capítulos. Los datos de la primera parte son simulados a propósito: así controlas la forma de la relación y puedes dibujar la curva verdadera debajo de las observaciones, sin descargar archivos externos.

# Bloque inicial del notebook: preparación del entorno
# En esta celda cargamos las librerías base de la serie y definimos una
# paleta visual consistente para que todas las gráficas se lean igual.
# Además, fijamos la semilla para que los datos sintéticos sean
# reproducibles en cada ejecución del notebook.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, SplineTransformer, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

SEMILLA = 42
rng = np.random.default_rng(SEMILLA)

plt.rcParams["figure.dpi"] = 110
plt.rcParams["axes.spines.top"] = False
plt.rcParams["axes.spines.right"] = False

COLOR_DATOS = "#4C78A8"
COLOR_VERDAD = "#333333"
COLOR_LINEAL = "#E45756"
COLOR_POLI = "#F58518"
COLOR_SPLINE = "#54A24B"
COLOR_KNOT = "#B279A2"
COLORES_BASE = [
    "#4C78A8",
    "#F58518",
    "#E45756",
    "#54A24B",
    "#B279A2",
    "#72B7B2",
    "#ECA400",
    "#9D755D",
]


def metricas_regresion(y_true, y_pred, nombre=""):
    mae = mean_absolute_error(y_true, y_pred)
    rmse = float(np.sqrt(mean_squared_error(y_true, y_pred)))
    r2 = r2_score(y_true, y_pred)
    encabezado = nombre if nombre else "Modelo"
    print(f"{encabezado}")
    print(f"  MAE  = {mae:.3f}")
    print(f"  RMSE = {rmse:.3f}")
    print(f"  R²   = {r2:.3f}")
    return mae, rmse, r2
# Muestra independiente para la selección y el reporte final.
rng_evaluacion = np.random.default_rng(SEMILLA + 1)


# Versiones del entorno usado para ejecutar este capítulo.
import sys
import importlib.metadata as metadatos
print('Python:', sys.version.split()[0])
for paquete in ['numpy', 'pandas', 'matplotlib', 'scipy', 'scikit-learn', 'statsmodels']:
    print(paquete + ':', metadatos.version(paquete))
import matplotlib.pyplot as plt
from cycler import cycler
plt.rcParams['axes.prop_cycle'] = cycler(color=['#1786B3', '#EE5400', '#487A35', '#7851A9']) + cycler(linestyle=['-', '--', '-.', ':'])
plt.rcParams['font.size'] = 10

Python: 3.13.14
numpy: 2.5.1
pandas: 3.0.5
matplotlib: 3.11.1
scipy: 1.18.0
scikit-learn: 1.9.0
statsmodels: 0.14.6

Las funciones de dibujo

El entorno queda con Python 3.13.14, NumPy 2.5.1, pandas 3.0.5 y scikit-learn 1.9.0.

Las funciones de la celda siguiente no son un alias de plot o scatter: cada una arma la figura completa (abre el lienzo, pinta lo que toca, pone ejes y la muestra), y en las celdas del análisis solo quedan los datos, el modelo y una llamada.

Función	Qué figura completa genera
figura_nube	Observaciones y, si quieres, la curva verdadera
figura_curvas	Nube, curva verdadera, uno o varios ajustes y knots
figura_residuos	Residuos contra \(X\), con la línea en cero
figura_paneles	Varios paneles: la misma nube y una predicción distinta en cada uno
figura_esquema_knots	Dibujo didáctico de cortes y regiones
figura_histograma_knots	Histograma de \(X\) y, debajo, dónde caen distintos juegos de knots
figura_cv	RMSE de validación cruzada contra la complejidad, con el valor elegido

Los nombres de ejes por defecto son los del tráfico; si cambias de fenómeno, pásalos tú.

def _cerrar(fig, titulo_superior=None, y_sup=1.03):
    if titulo_superior:
        fig.suptitle(titulo_superior, y=y_sup)
    plt.tight_layout()
    plt.show()


def _nube(ax, x, y, etiqueta="Observaciones", tamano=22):
    ax.scatter(
        x, y,
        c=COLOR_DATOS, s=tamano, alpha=0.7,
        edgecolors="white", linewidths=0.35, zorder=3, label=etiqueta,
    )


def _verdad(ax, x, y, etiqueta="Relación verdadera"):
    ax.plot(x, y, color=COLOR_VERDAD, lw=1.5, ls="--", label=etiqueta)


def _marca_knots(ax, knots):
    for xi in np.asarray(knots, dtype=float).ravel():
        ax.axvline(xi, color=COLOR_KNOT, ls=":", lw=1.4)


def _eje(ax, titulo, xlabel, ylabel=None, xlim=None):
    ax.set_title(titulo)
    ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if xlim is not None:
        ax.set_xlim(*xlim)
    ax.grid(True, alpha=0.25)


def _trazo(ax, x, y, color, etiqueta=None, trozos=None):
    if trozos is None:
        ax.plot(x, y, color=color, lw=2.2, label=etiqueta)
        return
    for i, (izq, der, *_) in enumerate(trozos):
        mascara = (x >= izq) & (x <= der)
        ax.plot(x[mascara], y[mascara], color=color, lw=2.2, label=etiqueta if i == 0 else None)


def figura_nube(
    x, y,
    x_linea=None, y_linea=None, etiqueta_linea="Relación verdadera",
    titulo="", xlabel="Hora del día", ylabel="Flujo (vehículos / min)", xlim=(0, 24),
):
    """Figura completa: la nube de observaciones y, si se pasan, una curva (casi siempre la verdadera)."""
    fig, ax = plt.subplots(figsize=(9.5, 4.8))
    _nube(ax, x, y, tamano=28)
    if x_linea is not None:
        _verdad(ax, x_linea, y_linea, etiqueta=etiqueta_linea)
    _eje(ax, titulo, xlabel, ylabel, xlim)
    ax.legend(frameon=False)
    _cerrar(fig)


def figura_curvas(
    x, y, x_grid, curvas,
    y_verdad=None, knots=None, trozos=None, colores=None,
    titulo="", xlabel="Hora del día", ylabel="Flujo (vehículos / min)", xlim=(0, 24),
    loc_leyenda="best",
):
    """Figura completa: observaciones + curva verdadera + uno o varios ajustes.

    curvas es un diccionario: {"Regresión lineal": prediccion_en_x_grid, ...}.
    knots dibuja líneas verticales punteadas. trozos parte el trazo si las piezas pueden saltar.
    """
    fig, ax = plt.subplots(figsize=(9.5, 4.8))
    _nube(ax, x, y)
    if y_verdad is not None:
        _verdad(ax, x_grid, y_verdad)
    if colores is None:
        colores = [COLOR_LINEAL, COLOR_POLI, COLOR_SPLINE, COLOR_KNOT] + COLORES_BASE
    for i, (nombre, pred) in enumerate(curvas.items()):
        _trazo(ax, x_grid, pred, colores[i % len(colores)], etiqueta=nombre, trozos=trozos)
    if knots is not None:
        _marca_knots(ax, knots)
    _eje(ax, titulo, xlabel, ylabel, xlim)
    ax.legend(frameon=False, loc=loc_leyenda)
    _cerrar(fig)


def figura_residuos(
    x, residuos,
    titulo="", xlabel="Hora del día", ylabel="Residuo (observado − predicción)", xlim=(0, 24),
):
    """Figura completa de residuos: puntos contra X y una línea horizontal en cero."""
    fig, ax = plt.subplots(figsize=(9.5, 4.2))
    _nube(ax, x, residuos, etiqueta=None, tamano=28)
    ax.axhline(0, color=COLOR_VERDAD, lw=1.2)
    _eje(ax, titulo, xlabel, ylabel, xlim)
    _cerrar(fig)


def figura_paneles(
    x, y, x_grid, paneles,
    y_verdad=None, knots=None, trozos=None, colores=None,
    titulo_figura=None, xlabel="Hora del día", ylabel="Flujo (vehículos / min)", xlim=(0, 24),
    n_cols=None, knots_por_panel=None,
):
    """Varios paneles con la misma nube. paneles = {"Grado 3": pred, "Grado 8": pred, ...}.

    En cada panel se pinta una predicción. knots_por_panel, si se pasa, es una lista de
    cortes distinta por panel (por ejemplo 3, 5, 8 y 15 knots).
    """
    nombres = list(paneles.keys())
    preds = list(paneles.values())
    n = len(nombres)
    if n_cols is None:
        n_cols = n if n <= 3 else 2
    n_filas = int(np.ceil(n / n_cols))
    fig, axes = plt.subplots(
        n_filas, n_cols,
        figsize=(4.2 * n_cols + 0.4, 3.7 * n_filas + 0.5),
        sharey=True, sharex=True,
    )
    axes_lista = np.atleast_1d(axes).ravel()
    if colores is None:
        colores = [COLOR_POLI] * n
    for i, ax in enumerate(axes_lista):
        if i >= n:
            ax.axis("off")
            continue
        _nube(ax, x, y, etiqueta=None, tamano=14)
        if y_verdad is not None:
            _verdad(ax, x_grid, y_verdad, etiqueta=None)
        _trazo(ax, x_grid, preds[i], colores[i % len(colores)], trozos=trozos)
        cortes = knots
        if knots_por_panel is not None:
            cortes = knots_por_panel[i]
        if cortes is not None:
            _marca_knots(ax, cortes)
        _eje(ax, nombres[i], xlabel, ylabel if i % n_cols == 0 else None, xlim)
    _cerrar(fig, titulo_superior=titulo_figura)


def figura_esquema_knots(x, y, knots, nombres_region, titulo=""):
    """Figura didáctica: una curva suave, líneas verticales en los knots y el nombre de cada región."""
    fig, ax = plt.subplots(figsize=(9.8, 4.6))
    ax.plot(x, y, color=COLOR_SPLINE, lw=2.6, zorder=3)
    knots = np.asarray(knots, dtype=float).ravel()
    for i, xi in enumerate(knots, start=1):
        ax.axvline(xi, color=COLOR_KNOT, ls=":", lw=1.8)
        ax.text(xi, y.max() + 3.5, rf"$\xi_{i}$", color=COLOR_KNOT, ha="center", fontsize=13, fontweight="bold")
    bordes = np.concatenate([[x.min()], knots, [x.max()]])
    for i, nombre in enumerate(nombres_region):
        ax.text(0.5 * (bordes[i] + bordes[i + 1]), y.min() - 4.5, nombre, ha="center", va="top", fontsize=11, color=COLOR_VERDAD)
    ax.set_xlabel("X")
    ax.set_ylabel("f(X)")
    ax.set_title(titulo)
    ax.set_xlim(0, 24)
    ax.set_ylim(y.min() - 10, y.max() + 10)
    ax.set_yticks([])
    ax.grid(True, alpha=0.25, axis="x")
    _cerrar(fig)


def figura_histograma_knots(
    x, grupos_knots, n_bins=24,
    titulo_hist="", titulo_knots="", xlabel="Hora del día", xlim=(0, 24),
):
    """Dos paneles: arriba un histograma de X; abajo los mismos puntos y varios juegos de knots.

    grupos_knots = {"Uniformes": (valores, color, estilo), "Cuantiles": (...)}.
    """
    fig, axes = plt.subplots(2, 1, figsize=(9.5, 6.2), sharex=True)
    axes[0].hist(x, bins=n_bins, color=COLOR_DATOS, edgecolor="white")
    axes[0].set_ylabel("Cantidad de mediciones")
    axes[0].set_title(titulo_hist)
    axes[0].grid(True, alpha=0.25, axis="y")
    axes[1].scatter(x, np.zeros_like(x), c=COLOR_DATOS, s=14, alpha=0.45, label="Observaciones (eje X)")
    for nombre, (xs, color, estilo) in grupos_knots.items():
        for i, xi in enumerate(np.asarray(xs, dtype=float).ravel()):
            axes[1].axvline(xi, color=color, ls=estilo, lw=1.6, label=nombre if i == 0 else None)
    axes[1].set_yticks([])
    _eje(axes[1], titulo_knots, xlabel, xlim=xlim)
    axes[1].legend(frameon=False, loc="upper left")
    _cerrar(fig)


def figura_cv(valores, rmse, std, elegido=None, xlabel="n_knots", ylabel="RMSE de validación cruzada", titulo=""):
    """Figura completa: error de CV contra la complejidad, con barras de 1 desv. estándar y el valor elegido."""
    fig, ax = plt.subplots(figsize=(9.5, 4.8))
    ax.errorbar(
        valores, rmse, yerr=std, fmt="o-",
        color=COLOR_SPLINE, ecolor=COLOR_KNOT, capsize=4, lw=2.0, ms=8,
        label="RMSE CV (media ± 1 desv. estándar)",
    )
    if elegido is not None:
        idx = list(valores).index(elegido)
        ax.scatter(
            [elegido], [rmse[idx]], s=160, facecolors="none",
            edgecolors=COLOR_LINEAL, lw=2.2, zorder=5,
            label=f"Seleccionado: {xlabel} = {elegido}",
        )
        ax.axvline(elegido, color=COLOR_LINEAL, ls="--", lw=1.2, alpha=0.85)
        ax.annotate(
            "Mejor complejidad según CV",
            xy=(elegido, rmse[idx]),
            xytext=(elegido - 3.2, rmse[idx] + 3.5),
            fontsize=10, color=COLOR_LINEAL,
            arrowprops=dict(arrowstyle="->", color=COLOR_LINEAL, lw=1.2),
        )
    _eje(ax, titulo, xlabel, ylabel)
    ax.legend(frameon=False)
    _cerrar(fig)

20.2 Cuando una recta no alcanza

Construyes un día de tráfico de 24 horas con 180 mediciones. La relación verdadera tiene dos picos, la entrada y la salida laborales, sobre un ritmo de fondo que sube y baja con suavidad, y encima hay ruido: no todos los martes se parecen, hay un accidente, un partido, una lluvia.

Tres cosas importan de este diseño. La curva no es una recta ni un arco suave de un solo máximo. Conoces la función que genera los datos, así que más adelante podrás juzgar si el modelo recuperó la forma o solo recorrió los puntos. Y el ruido existe, de modo que un modelo que pase exactamente por cada observación estará memorizando el martes, no el patrón.

La pregunta de esta sección es todavía más simple que elegir un modelo: si una recta puede contar esta historia.

Generar el día de tráfico y mirar los puntos

Generas 180 mediciones a lo largo del día y dibujas la curva verdadera solo como referencia didáctica; en un problema real esa curva no la tienes, y por eso el modelo tiene que inferirla.

Las primeras figuras y métricas son ilustrativas e in-sample, sobre esta muestra de 180 mediciones. En la misma celda generas una segunda muestra independiente y reservas su prueba antes de inspeccionarla, y solo esa segunda muestra interviene en la selección por validación cruzada y en el reporte final.

# Generamos un dataset sintético para representar un patrón de tráfico con dos horas pico.
# La idea es construir una relación no lineal y luego añadir ruido para simular variabilidad real.
# La función 'flujo_verdadero' representa la curva teórica del fenómeno, pero en la práctica
# no la vemos: el modelo debe inferirla a partir de las observaciones ruidosas.

def flujo_verdadero(hora):
    """Vehículos por minuto en función de la hora del día."""
    # Componente suave de fondo: el flujo sube y baja durante el día.
    ritmo_base = 24 + 7 * np.sin((hora - 5) * np.pi / 12)
    # Dos picos locales: uno por la mañana y otro por la tarde.
    pico_manana = 80 * np.exp(-0.5 * ((hora - 8.0) / 0.95) ** 2)
    pico_tarde = 92 * np.exp(-0.5 * ((hora - 18.2) / 1.15) ** 2)
    return ritmo_base + pico_manana + pico_tarde


# 180 observaciones distribuidas uniformemente en el rango del día, con ruido gaussiano.
# La variable 'hora' será la entrada X, y 'flujo' la variable objetivo y.
n = 180
hora = np.sort(rng.uniform(0.4, 23.6, size=n))
flujo = flujo_verdadero(hora) + rng.normal(0, 6.5, size=n)

# Preparamos las matrices para trabajar con sklearn.
# En este caso, X tiene una sola columna: la hora del día.
X = hora.reshape(-1, 1)
y = flujo.copy()

# Una segunda muestra del mismo mecanismo queda separada del material ilustrativo.
hora_eval = np.sort(rng_evaluacion.uniform(0.4, 23.6, size=n))
X_eval = hora_eval.reshape(-1, 1)
y_eval = flujo_verdadero(hora_eval) + rng_evaluacion.normal(0, 6.5, size=n)
X_ent, X_prueba, y_ent, y_prueba = train_test_split(
    X_eval, y_eval, test_size=0.25, random_state=SEMILLA
)

# Creamos una grilla fina para dibujar la curva de predicción sobre todo el eje X.
x_grid = np.linspace(hora.min(), hora.max(), 400)
X_grid = x_grid.reshape(-1, 1)
y_real_grid = flujo_verdadero(x_grid)

# Resumen del dataset generado.
print(f"Observaciones: {n}")
print(f"Hora mínima: {hora.min():.2f}  |  hora máxima: {hora.max():.2f}")
print(f"Flujo medio: {y.mean():.1f} vehículos/min")

figura_nube(
    hora, y,
    x_linea=x_grid, y_linea=y_real_grid,
    titulo="El flujo vehicular cambia de régimen a lo largo del día",
)

Observaciones: 180
Hora mínima: 0.57  |  hora máxima: 23.42
Flujo medio: 46.0 vehículos/min

La figura de las 180 mediciones ilustrativas muestra dos picos, cerca de las 8:00 y las 18:12, separados por un valle. Si tu modelo solo puede dibujar una recta, va a tener que atravesar las dos montañas como si fueran ruido alrededor de una media. El patrón está a la vista; la recta no tiene permiso para seguirlo.

Esa es la distinción que importa: los datos no son “ruidosos y sin estructura”. Tienen una estructura que el modelo lineal no puede nombrar.

Ajustar la recta que ya conoces

Usas la misma herramienta del capítulo de regresión lineal del nivel Básico, \(y = \beta_0 + \beta_1 x + \varepsilon\), con las métricas de siempre. Los residuos frente a la hora te van a mostrar en qué tramos se concentra el error de la recta.

modelo_lineal = LinearRegression()  # instancia el modelo
modelo_lineal.fit(X, y)
pred_lineal = modelo_lineal.predict(X)
pred_lineal_grid = modelo_lineal.predict(X_grid)

print(
    f"Intercepto = {modelo_lineal.intercept_:.2f}  |  "
    f"pendiente = {modelo_lineal.coef_[0]:.3f} vehículos/min por hora"
)
metricas_regresion(y, pred_lineal, "Regresión lineal (ajuste in-sample)")

figura_curvas(
    hora, y, x_grid,
    curvas={"Regresión lineal": pred_lineal_grid},
    y_verdad=y_real_grid,
    colores=[COLOR_LINEAL],
    titulo="Una recta no distingue las horas pico",
)

Intercepto = 31.89  |  pendiente = 1.192 vehículos/min por hora
Regresión lineal (ajuste in-sample)
  MAE  = 24.500
  RMSE = 30.217
  R²   = 0.062

Sobre las 180 mediciones ilustrativas, la recta ajustada tiene la forma

\[ \hat{y} = \beta_0 + \beta_1 x \]

Los números concretos salen en la celda anterior: un intercepto de 31.89 vehículos por minuto y una pendiente de 1.192 vehículos por minuto por hora. Sobre el rango observado (de las 0:34 a las 23:25), esa pendiente hace subir la recta unos 27 vehículos por minuto, de ~33 a ~60. Eso resume todo el día en una sola tendencia global, que sube porque el pico de la tarde es más alto y más ancho que el de la mañana. El flujo sube hacia cada pico y vuelve a bajar después; esas subidas y bajadas se compensan en buena parte, y la pendiente solo recoge esa tendencia leve. La pendiente no es incorrecta de forma matemática; lo que está mal es la forma funcional. El modelo no está “detectando que la hora casi no importa”: está obligado a decir que importa siempre igual, y con esa obligación lo único que puede contar es una tendencia promedio. La relación real tiene un fondo suave, un pico de la mañana y un pico de la tarde, y una sola recta no puede representar esos cambios locales sin promediarlos.

Las métricas que imprimiste cuentan la misma historia:

• El MAE, cerca de 25, dice que en promedio cada predicción se equivoca en unos 25 vehículos por minuto. Como la variable objetivo está en esas mismas unidades, el número se lee directo. Un error típico de ese tamaño significa que, incluso en el mejor caso del ajuste global, la recta falla en una fracción grande del flujo de cada instante.
• El RMSE, cerca de 30, refleja que el error cuadrático amplifica los errores más grandes, por eso nunca es menor que el MAE. Esos errores grandes aparecen en los picos de la mañana y la tarde, donde la recta subestima o sobreestima mucho más que el promedio. El RMSE muestra el peso de los errores grandes, que aquí se concentran en los picos.
• El \(R^2\), cerca de 0.06, indica que la recta solo explica alrededor del 6% de la variación total. Es decir, la recta no captura la estructura del fenómeno y predice casi igual que la media del flujo (una línea horizontal). Un valor tan bajo, por sí solo, no dice por qué: puede deberse a mucho ruido o a una forma equivocada. Aquí el gráfico ya mostró que es lo segundo, y los residuos lo van a confirmar.

En pocas palabras: MAE y RMSE miden cuánto se equivoca el modelo; \(R^2\) mide cuánto de la variación real logra explicar. Aquí se equivoca bastante y explica muy poco. El \(R^2\) de 0.062 refleja aquí una forma que representa mal los dos picos simulados; el gráfico te dice dónde se concentra el error: en los dos picos, exactamente donde a un operador de tráfico le importaría acertar.

En los capítulos de evaluación y de validación cruzada aprendiste que un error grande puede ser sobreajuste, subajuste o simple azar de la partición. Aquí todavía ni siquiera saliste del ajuste dentro de la muestra (in-sample: el error se mide sobre los mismos datos con los que se ajustó) y el error ya es enorme. La recta no falla porque el ruido sea enorme ni porque hayas evaluado mal: falla porque no es el objeto correcto.

Los residuos cuentan la misma historia

Si el modelo estuviera bien especificado, los residuos no deberían guardar la forma de las horas pico. La figura siguiente lo comprueba.

residuos_lineal = y - pred_lineal  # residuo = observado − predicción

figura_residuos(
    hora, residuos_lineal,
    titulo="Los residuos de la recta todavía tienen las dos horas pico",
)

Los residuos no forman una nube aleatoria alrededor de cero. Tienen una estructura clara en función de la hora: son positivos en los dos picos, negativos en la madrugada, el mediodía y la noche, y el signo cambia justo donde cambia el patrón del tráfico. La recta no está ignorando el tráfico: está forzando una forma demasiado rígida y deja toda la estructura en los residuos, que repiten la forma de la señal que quieres explicar. Ese patrón coincide con la forma que la recta no representa y motiva ampliar la función. En las siguientes figuras evaluarás formas que pueden representar la curvatura de esos dos picos.

Un buen modelo deja residuos sin patrón: dispersos alrededor de cero, sin tendencia, sin curvatura ni picos repetidos. Si los residuos siguen la forma del fenómeno original, entonces el modelo está dejando la estructura en el error y no está capturando bien la relación.

20.3 Regresión polinómica: curvar sin salir de la regresión lineal

La primera ampliación natural es permitir potencias de \(X\):

\[ y = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots + \beta_p x^p + \varepsilon \]

Fíjate en lo que no cambió. Los parámetros \(\beta_0,\ldots,\beta_p\) siguen entrando de forma lineal. Si llamas \(z_1 = x\), \(z_2 = x^2\), \(\ldots\), \(z_p = x^p\), esto es una regresión lineal múltiple sobre variables derivadas (features). Por eso se estima con el mismo LinearRegression de siempre: el truco está en expandir \(X\), no en cambiar el motor de estimación.

Lo que sí cambió es la forma que el modelo puede dibujar. Grado 1 es la recta. Grado 2 puede hacer un arco. Grado 3 puede tener un máximo y un mínimo. Cada grado extra permite un cambio de dirección más.

Para comparar las curvas, pregunta: ¿qué está haciendo realmente este modelo en todo el dominio a la vez?

Grado	Forma típica	Qué puede representar
1	línea recta	tendencia global constante
2	arco	una sola curvatura, como una U o una montaña
3	dos cambios de dirección	un máximo y un mínimo
4	hasta tres cambios de dirección	por ejemplo, dos máximos con un mínimo entre ellos
5	hasta cuatro cambios de dirección	más detalle, y también más riesgo de oscilar donde no debe

Este cuadro resume la idea: a medida que aumenta el grado, la curva gana libertad para seguir la forma real de los datos, pero también se vuelve más sensible al ruido.

Mira un ejemplo pequeño de esa transformación. Si empiezas con una sola variable \(x\), PolynomialFeatures crea nuevas columnas con \(x\), \(x^2\) y \(x^3\) para que el modelo pueda aprender una curva.

El parámetro degree dice hasta qué potencia se expande la variable: con degree=2 la base es \([x, x^2]\), y con degree=3, \([x, x^2, x^3]\). Un grado mayor da más libertad para doblar la curva y también aumenta el riesgo de sobreajustar. include_bias=False evita añadir una columna de unos, porque el modelo ya incluye su propio intercepto y la duplicaría.

from sklearn.preprocessing import PolynomialFeatures

x = np.array([1, 2, 3, 4]).reshape(-1, 1)

for grado in [2, 3]:
    poly = PolynomialFeatures(degree=grado, include_bias=False)
    X_poly = poly.fit_transform(x)
    print(f"Grado {grado}:")
    print(np.round(X_poly, 3))
    print()

# Con grado 2 cada fila se vuelve [x, x^2].
# Con grado 3 cada fila se vuelve [x, x^2, x^3].
# Más potencia = más flexibilidad para capturar curvatura, pero también más riesgo de
# ajustar ruido en lugar de tendencia real.

Grado 2:
[[ 1.  1.]
 [ 2.  4.]
 [ 3.  9.]
 [ 4. 16.]]

Grado 3:
[[ 1.  1.  1.]
 [ 2.  4.  8.]
 [ 3.  9. 27.]
 [ 4. 16. 64.]]

20.4 El problema de una única función para todo el dominio

Un polinomio global usa la misma fórmula para todo el dominio de \(x\). Eso es importante porque si la relación real tiene dos picos distintos, ese mismo modelo debe intentar acomodarlos con la misma forma matemática en todas partes.

A medida que sube el grado, la curva gana flexibilidad, pero también puede volverse más inestable. No es una falla de programación ni un error de Python: es la naturaleza de una sola función global tratando de explicar varias zonas distintas del dominio.

Subir el grado con calma: 1, 2 y 3

Tres paneles, no todas las curvas encimadas en uno. En cada uno, la misma nube, la misma curva verdadera y un solo polinomio. PolynomialFeatures fabrica las columnas \(x, x^2, \ldots\); LinearRegression estima los \(\beta\). El StandardScaler, que va primero, no cambia la idea: reduce los problemas numéricos de elevar a potencias altas una hora que va de 0 a 24.

El Pipeline de la celda siguiente encadena los tres pasos en un estimador. Escala, expande a potencias y estima los \(\beta\) como un solo objeto. Así, fit y predict recorren la misma cadena y no se te olvida escalar a la hora de predecir. La regresión sigue siendo lineal en los coeficientes y ahora actúa sobre la expansión polinómica de la entrada. Eso es lo que permite doblar la curva sin cambiar el algoritmo base.

Si la entrada original es una sola columna \(x\), entonces para grado \(p\) la matriz construida es:

\[ X_p = \begin{bmatrix} x_1 & x_1^2 & x_1^3 & \cdots & x_1^p \\ x_2 & x_2^2 & x_2^3 & \cdots & x_2^p \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_n & x_n^2 & x_n^3 & \cdots & x_n^p \end{bmatrix} \]

Y la predicción sigue siendo una regresión lineal sobre esa matriz:

\[ \hat{y} = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots + \beta_p x^p \]

La representación de la entrada cambia al incorporar columnas de potencias. Antes tenías una sola columna; ahora tienes varias que codifican la curvatura. Esa es la idea que permite aproximar una relación no lineal con una estructura lineal en los parámetros.

# Tres grados, uno por uno. No buscamos un ganador: queremos ver qué forma
# puede cada uno y si alcanza para dos picos. Los tres son UNA fórmula
# para todo el día.

# El Pipeline es una tubería: cada paso recibe la salida del anterior.
#
#   hora (X)
#      → StandardScaler          (x − media) / desv.   →  X_escalada
#      → PolynomialFeatures      [x, x², …, xᵖ]        →  X_polinómica
#      → LinearRegression        estima los β          →  ŷ
#
# PolynomialFeatures no ve la hora cruda (0 a 24): ve la hora ya escalada.
# Si x³ se fabricara antes, 24³ = 13824: columnas de tamaños tan distintos
# vuelven inestable el cálculo.
# include_bias=False: LinearRegression ya trae el intercepto.
#
# fit recorre la cadena y cada paso aprende lo suyo (media y desv.; luego β).
# predict recorre la misma cadena con datos nuevos.


paneles = {}  # título de cada panel → su predicción en la malla


# Grado 1: una recta. Una pendiente para las 24 horas.
modelo_g1 = Pipeline(
    [
        ("escala", StandardScaler()),  # centra y escala la hora
        ("poly", PolynomialFeatures(degree=1, include_bias=False)),  # columnas: x
        ("modelo", LinearRegression()),  # estima β0 y β1
    ]
)
modelo_g1.fit(X, y)
_, rmse_g1, r2_g1 = metricas_regresion(
    y, modelo_g1.predict(X), "Polinomio grado 1 (in-sample)"
)
paneles[f"Grado 1  |  RMSE={rmse_g1:.1f}  $R^2$={r2_g1:.2f}"] = modelo_g1.predict(X_grid)


# Grado 2: un arco. Tiene a lo sumo un pico, así que no puede formar los dos.
modelo_g2 = Pipeline(
    [
        ("escala", StandardScaler()),  # misma hora escalada que en el grado 1
        ("poly", PolynomialFeatures(degree=2, include_bias=False)),  # columnas: x, x²
        ("modelo", LinearRegression()),  # estima β0, β1 y β2
    ]
)
modelo_g2.fit(X, y)
_, rmse_g2, r2_g2 = metricas_regresion(
    y, modelo_g2.predict(X), "Polinomio grado 2 (in-sample)"
)
paneles[f"Grado 2  |  RMSE={rmse_g2:.1f}  $R^2$={r2_g2:.2f}"] = modelo_g2.predict(X_grid)


# Grado 3: en este rango, un solo máximo. Se acerca, pero sigue siendo un cubo
# global: no trata cada hora pico por separado.
modelo_g3 = Pipeline(
    [
        ("escala", StandardScaler()),  # otra vez: escala primero, potencias después
        ("poly", PolynomialFeatures(degree=3, include_bias=False)),  # columnas: x, x², x³
        ("modelo", LinearRegression()),  # estima β0, β1, β2 y β3
    ]
)
modelo_g3.fit(X, y)
_, rmse_g3, r2_g3 = metricas_regresion(
    y, modelo_g3.predict(X), "Polinomio grado 3 (in-sample)"
)
paneles[f"Grado 3  |  RMSE={rmse_g3:.1f}  $R^2$={r2_g3:.2f}"] = modelo_g3.predict(X_grid)


figura_paneles(
    hora,  # eje X de las observaciones (la hora)
    y,  # eje Y de las observaciones (el flujo)
    x_grid,  # grilla fina: ahí se dibujan las curvas, no en los puntos
    paneles,  # título de cada cuadro y su predicción ya calculada
    y_verdad=y_real_grid,  # relación verdadera, punteada, para comparar formas
    colores=[COLOR_POLI] * 3,  # el mismo color en los tres polinomios
    titulo_figura="Polinomios de grado bajo: todavía una sola onda para dos picos",  # título general
)

Polinomio grado 1 (in-sample)
  MAE  = 24.500
  RMSE = 30.217
  R²   = 0.062
Polinomio grado 2 (in-sample)
  MAE  = 23.349
  RMSE = 28.663
  R²   = 0.156
Polinomio grado 3 (in-sample)
  MAE  = 22.881
  RMSE = 28.418
  R²   = 0.170

El grado 1 es la recta de antes (RMSE 30.2, \(R^2\) = 0.06). El grado 2 añade un arco, y con un solo máximo no puede servir a los dos picos: lo coloca hacia las 14 h, justo en el valle entre ellos (RMSE 28.7, \(R^2\) = 0.16). El grado 3 podría subir, bajar y volver a subir, pero en estos datos no le conviene: su curva tiene un único máximo, hacia las 15.5 h, y apenas mejora (RMSE 28.4, \(R^2\) = 0.17). Sigue siendo una función polinómica global, con muy pocos cambios de curvatura permitidos. Los dos picos afilados del tráfico no son un cubo suave. El modelo está intentando contar dos historias locales con un vocabulario demasiado pequeño.

Qué ocurre si aumentas el grado

Grados 5, 10 y 12. Más parámetros, más permiso para doblarse. La tentación es creer que “más grado = más verdad”. Mira los extremos del día, no solo el centro.

# En la celda anterior escribiste la tubería a mano (grado 1, 2 y 3) para
# ver cada paso. A partir de aquí la misma tubería se usa muchas veces:
# grados altos, k-fold y comparación de modelos.
# Por claridad y reutilización, la metemos en una función: le pasas el
# grado y te devuelve el Pipeline listo para .fit() / .predict().
def pipeline_polinomico(grado):
    return Pipeline(
        [
            ("escala", StandardScaler()),
            ("poly", PolynomialFeatures(degree=grado, include_bias=False)),
            ("modelo", LinearRegression()),
        ]
    )


# Tres grados altos. El for solo cambia el número; el resto es idéntico.
grados_altos = [5, 10, 12]
paneles = {}
for grado in grados_altos:
    modelo = pipeline_polinomico(grado)
    modelo.fit(X, y)
    _, rmse, r2 = metricas_regresion(y, modelo.predict(X), f"Polinomio grado {grado} (in-sample)")
    paneles[f"Grado {grado}  |  RMSE={rmse:.1f}  $R^2$={r2:.2f}"] = modelo.predict(X_grid)

figura_paneles(
    hora, y, x_grid, paneles,
    y_verdad=y_real_grid,
    colores=[COLOR_POLI] * 3,
    titulo_figura="Más grado: mejor centro, peor conducta en los extremos",
)

Polinomio grado 5 (in-sample)
  MAE  = 19.006
  RMSE = 23.779
  R²   = 0.419
Polinomio grado 10 (in-sample)
  MAE  = 11.269
  RMSE = 13.663
  R²   = 0.808
Polinomio grado 12 (in-sample)
  MAE  = 9.619
  RMSE = 11.645
  R²   = 0.861

Al subir el grado, el polinomio gana flexibilidad y el error de entrenamiento baja: RMSE in-sample de 23.8 (grado 5), 13.7 (grado 10) y 11.6 (grado 12). Eso no es gratis. Un polinomio de grado alto es una función global: los mismos \(\beta\) tienen que servir para las 5 de la mañana y para las 6 de la tarde. Para clavar un pico afilado, la función necesita ondular; esas ondulaciones no se quedan en el pico, viajan por el resto del dominio. El síntoma típico es el vaivén en los extremos, donde la curva solo tiene vecinos de un lado que la anclen: el grado 5 se desploma hasta −45.6 vehículos por minuto a las 23:25, un flujo negativo, imposible, y el grado 12 arranca en 38 a las 0:34, cuando la curva verdadera está cerca de 18.

Estás viendo tres regímenes de complejidad, los mismos de los capítulos de regularización y de evaluación:

• subajuste: grados 1–3, la forma no alcanza;
• ajuste razonable: un grado intermedio se acerca, con concesiones;
• exceso de flexibilidad: la curva se vuelve inestable, sobre todo en los extremos.

Conviene separar tres cosas que se mezclan con facilidad: que la forma no alcance, que la curva oscile en los bordes y que el modelo ajuste el ruido. La figura muestra las dos primeras: el grado 12 todavía no alcanza los picos (su RMSE in-sample, 11.6, sigue por encima del ruido con que se generaron los datos, 6.5) y ya oscila en los extremos. La tercera puede estar ocurriendo a la vez; con el ajuste in-sample no se puede afirmar ni descartar.

Los RMSE in-sample de 23.8, 13.7 y 11.6 deben contrastarse mediante validación: reducirlos puede introducir errores en otras regiones del eje \(X\).

20.5 El fenómeno de Runge: el polinomio oscila donde la función no lo hace

Eso que acabas de ver en el tráfico no es una falla de Python ni un accidente de estos datos. Está emparentado con un resultado clásico: el fenómeno de Runge.

Carl Runge mostró, a comienzos del siglo XX, un hecho incómodo: si interpolas ciertas funciones suaves con un polinomio de grado cada vez más alto, sobre puntos igualmente espaciados, no es cierto que te acerques a la función. En el centro el ajuste puede mejorar, mientras que cerca de los extremos la curva oscila y, al subir el grado, esas oscilaciones crecen. Más grado no implica más verdad: puede implicar peor conducta justo en los bordes. No ocurre con todas las funciones ni con cualquier forma de elegir los puntos, pero sí con funciones tan inocentes como la del ejemplo de abajo. En el tráfico la situación es distinta: no se interpola, hay ruido, los puntos no están equiespaciados y hay muchos más datos que parámetros. Tómalo como una analogía que ayuda a entender las oscilaciones de borde, no como el diagnóstico de lo que pasó allí.

La razón es geométrica, no computacional. Un polinomio de grado \(n\) es una sola función global: los mismos coeficientes \(\beta_0,\ldots,\beta_n\) tienen que servir en todo el intervalo. Para pasar cerca de muchos puntos del centro, la curva se ve obligada a compensar en los extremos, donde hay menos vecinos que la anclen. Esas ondulaciones no las pidió \(f\): las inventa el modelo.

Un ejemplo clásico aísla la idea sin el ruido del tráfico. La función verdadera es suave: no tiene saltos, ni picos afilados, ni oscilaciones propias. Es una campana centrada en 0, con máximo 1, que decae de forma monótona hacia los extremos (en \(x=\pm 1\) vale \(1/26 \approx 0.04\)).

\[ f(x) = \frac{1}{1 + 25x^2}, \qquad x \in [-1, 1] \]

En la celda siguiente tomas trece puntos igualmente espaciados de esta función, sin ruido, y les ajustas polinomios de grado 3, 8 y 12. El centro se ajusta bien; los bordes oscilan, y como la función es suave, cualquier oscilación en los bordes la inventó el modelo. El grado no es un volumen que se sube sin consecuencias.

x_runge = np.linspace(-1, 1, 13)
y_runge = 1.0 / (1.0 + 25.0 * x_runge**2)
x_runge_fino = np.linspace(-1, 1, 500)
y_runge_real = 1.0 / (1.0 + 25.0 * x_runge_fino**2)
X_runge = x_runge.reshape(-1, 1)
X_runge_fino = x_runge_fino.reshape(-1, 1)

paneles = {}
for grado in [3, 8, 12]:
    modelo = Pipeline(
        [
            ("poly", PolynomialFeatures(degree=grado, include_bias=False)),
            ("modelo", LinearRegression()),
        ]
    )
    modelo.fit(X_runge, y_runge)
    paneles[f"Grado {grado}"] = modelo.predict(X_runge_fino)

figura_paneles(
    x_runge, y_runge, x_runge_fino, paneles,
    y_verdad=y_runge_real,
    colores=[COLOR_POLI] * 3,
    titulo_figura="Fenómeno de Runge: el polinomio oscila donde la función no lo hace",
    xlabel="x",
    ylabel="f(x)",
    xlim=None,
)

El grado 3 se queda corto: aplana la campana y deja el pico en 0.49 cuando la función vale 1. El grado 8 se acerca en el centro (máximo 0.86) y ya muestra ondulaciones cerca de los bordes, que bajan hasta −0.09 aunque \(f\) nunca es negativa. El grado 12, con trece puntos, tiene exactamente un parámetro por observación (trece coeficientes, contando el intercepto): interpola los trece puntos sin error y, entre los dos últimos de cada lado, se hunde hasta −3.6, cuando la función verdadera no baja de 0.04. La oscilación entre observaciones se debe al interpolante y no aparece en la función generadora.

El ejemplo de Runge muestra un límite concreto: una sola función polinómica para todo el dominio acopla regiones que no deberían estar acopladas. El grado 12 interpola los trece puntos y aun así introduce errores grandes entre ellos.

Entonces la pregunta natural es otra.

20.6 Una idea más natural: dividir el problema

¿Y si, en lugar de obligar a una única función a explicar las 3 de la mañana y las 6 de la tarde, permites que diferentes regiones de \(X\) tengan comportamientos diferentes?

Eso es un polinomio por tramos: el dominio se corta en intervalos y en cada intervalo rige un polinomio propio. La flexibilidad deja de ser “subir el grado de todo el día” y pasa a ser “dónde permito que la historia cambie”.

Antes de ajustar nada, hay que nombrar esos cortes.

Qué es un knot

Los puntos de corte se llaman knots (nodos, puntos de corte). Se escriben

\[ \xi_1,\ \xi_2,\ \xi_3 \]

Con tres knots el eje \(X\) queda partido en cuatro regiones. La figura de abajo no busca ajustar tráfico: solo enseña esa geometría.

# Tres cortes → cuatro regiones. Esto es solo geometría: no ajustamos datos.
knots_esquema = np.array([6.0, 12.0, 18.0])

x_esquema = np.linspace(0, 24, 400)
# Curva ilustrativa, suave, sin saltos. Los knots no la "rompen": solo la marcan.
y_esquema = 40 + 18 * np.sin((x_esquema - 3) * np.pi / 12)

nombres_region = ["Región 1", "Región 2", "Región 3", "Región 4"]

figura_esquema_knots(
    x_esquema, y_esquema, knots_esquema,
    nombres_region=nombres_region,
    titulo="Un knot es un valor de X donde permitimos que cambie la forma local",
)

Los knots son valores de \(X\) donde permites que cambie la forma local de la función. Tres precisiones que no conviene olvidar:

• Un knot no tiene que coincidir con una observación. Es una decisión de diseño sobre el eje \(X\), no un punto especial de la muestra.
• Un knot no significa que la curva sea discontinua. Solo indica dónde el modelo puede cambiar de pieza.
• A la izquierda de \(\xi_k\) rige una forma; a la derecha, otra. Si esas formas se unen bien o mal es otro problema: el que viene ahora.

20.7 Polinomios por tramos

Con dos knots el modelo se escribe así:

\[ f(x)= \begin{cases} p_1(x), & x < \xi_1 \\[4pt] p_2(x), & \xi_1 \leq x < \xi_2 \\[4pt] p_3(x), & x \geq \xi_2 \end{cases} \]

Cada \(p_i(x)\) es un polinomio y gobierna solo su región. Los knots \(\xi_i\) separan esas regiones.

Si estimas cada \(p_i\) de manera independiente, nada garantiza que las piezas coincidan en los knots. Estimar piezas independientes permite flexibilidad local sin imponer continuidad en los cortes.

Volviendo al tráfico: los tres cortes del esquema eran solo para ver la geometría; para el experimento bastan dos knots, \(\xi_1=8\) y \(\xi_2=17\), que parten el día en madrugada, jornada y noche. En horas del día queda así:

\[ f(x)= \begin{cases} p_1(x), & x < 8 \quad \text{(madrugada)} \\[4pt] p_2(x), & 8 \leq x < 17 \quad \text{(jornada)} \\[4pt] p_3(x), & x \geq 17 \quad \text{(noche)} \end{cases} \]

No hay un modelo nuevo: son tres LinearRegression, cada una solo con los puntos de su intervalo. Después las pegas y miras el empalme.

knots_demo = np.array([8.0, 17.0])  # cortes: madrugada | jornada | noche

# Región 1: misma LinearRegression de siempre, solo con horas < 8
m1 = hora < 8  # máscara: True/False por observación
modelo_madrugada = LinearRegression()  # recta: y = β0 + β1 x
modelo_madrugada.fit(X[m1], y[m1])  # X[m1] se queda solo con las filas True

# Región 2: 8 ≤ hora < 17. Los paréntesis son obligatorios.
m2 = (hora >= 8) & (hora < 17)  # & = "y" casilla por casilla
modelo_jornada = LinearRegression()
modelo_jornada.fit(X[m2], y[m2])  # esta recta no ve madrugada ni noche

# Región 3: hora ≥ 17
m3 = hora >= 17
modelo_noche = LinearRegression()
modelo_noche.fit(X[m3], y[m3])

# Hasta aquí solo hay tres modelos ajustados; falta predecir en la malla para dibujarlos.

# Predicción en la grilla fina: cada hora usa SOLO la recta de su región
pred_tramos_grid = np.full(x_grid.shape, np.nan)  # vector vacío (NaN) del tamaño de la grilla
pred_tramos_grid[x_grid < 8] = modelo_madrugada.predict(X_grid[x_grid < 8])  # tramo 1
pred_tramos_grid[(x_grid >= 8) & (x_grid < 17)] = modelo_jornada.predict(
    X_grid[(x_grid >= 8) & (x_grid < 17)]
)  # tramo 2
pred_tramos_grid[x_grid >= 17] = modelo_noche.predict(X_grid[x_grid >= 17])  # tramo 3

# Lo mismo sobre las observaciones, para MAE / RMSE / R²
pred_tramos_obs = np.full(hora.shape, np.nan)
pred_tramos_obs[m1] = modelo_madrugada.predict(X[m1])
pred_tramos_obs[m2] = modelo_jornada.predict(X[m2])
pred_tramos_obs[m3] = modelo_noche.predict(X[m3])

piezas_lineales = [(0.0, 8.0), (8.0, 17.0), (17.0, 24.0)]  # para no unir el salto con una línea

figura_curvas(
    hora,  # eje X de las observaciones
    y,  # flujo observado
    x_grid,  # grilla donde se dibujan las tres rectas
    curvas={"Recta por tramos": pred_tramos_grid},  # las tres predicciones ya pegadas
    knots=knots_demo,  # verticales en 8 y 17
    trozos=piezas_lineales,  # pinta cada tramo aparte: así se ve el salto
    colores=[COLOR_POLI],
    titulo="Tres rectas independientes: la predicción puede saltar en el knot",
)

# Primero TODO el día: las tres predicciones pegadas, no una zona sola.
_ = metricas_regresion(y, pred_tramos_obs, "Tres rectas pegadas — todo el día (in-sample)")

# Ahora POR ZONA: cada recta solo contra los puntos de su intervalo.
_ = metricas_regresion(y[m1], pred_tramos_obs[m1], "Solo madrugada (hora < 8)")
_ = metricas_regresion(y[m2], pred_tramos_obs[m2], "Solo jornada (8 ≤ hora < 17)")
_ = metricas_regresion(y[m3], pred_tramos_obs[m3], "Solo noche (hora ≥ 17)")

Tres rectas pegadas — todo el día (in-sample)
  MAE  = 14.779
  RMSE = 17.556
  R²   = 0.683
Solo madrugada (hora < 8)
  MAE  = 14.700
  RMSE = 16.949
  R²   = 0.633
Solo jornada (8 ≤ hora < 17)
  MAE  = 14.564
  RMSE = 18.141
  R²   = 0.136
Solo noche (hora ≥ 17)
  MAE  = 15.190
  RMSE = 17.402
  R²   = 0.814

Revisar los empalmes

Cada modelo local puede ajustar bien su región, y la pregunta pasa a ser qué ocurre en los puntos donde una región termina y comienza la siguiente. En el gráfico, busca \(\xi_1\) y \(\xi_2\) y mira si la predicción llega al corte con el mismo valor desde la izquierda y desde la derecha: si un vehículo pasa de las 7:59 a las 8:01, el modelo no debería cambiarle el flujo de golpe.

Cada pieza puede ser razonable dentro de su región. El problema aparece en el empalme. En la figura, a las 8:00 la recta de la madrugada llega a ~76 vehículos por minuto y la de la jornada arranca en ~55; a las 17:00 el salto es mayor, de ~29 a ~120. Un salto significa que el modelo trata el knot como una pared: a las 7:59 una historia, a las 8:01 otra. Las métricas también muestran el costo de cortar mal: la recta de la jornada tiene \(R^2\) = 0.136, porque el knot de las 17 cae antes del pico de la tarde y la subida de las 16–17 h queda dentro de un tramo que la recta no puede seguir. En un fenómeno continuo eso no es flexibilidad, es un artefacto de estimar las piezas por separado.

La idea de partir el dominio era correcta. Faltaba una restricción: las piezas tienen que reunirse bien.

20.8 Splines: polinomios por tramos que se conectan suavemente

Un spline es un polinomio por tramos con condiciones de continuidad en los knots: las piezas no se estiman sueltas, están obligadas a coincidir en valor, pendiente y, si el grado alcanza, curvatura. El más usado, y el foco de este capítulo, es el cúbico, porque cada región es un polinomio de grado 3: alcanza para imponer continuidad de \(f\), \(f'\) y \(f''\) y la curva se siente de una pieza, sin las ondulaciones de un grado más alto.

Eso tiene nombre. Se dice que una función es \(C^k\) si es continua y se puede derivar \(k\) veces, y esas \(k\) derivadas también son continuas. El superíndice indica cuántas derivadas continuas se exigen en el knot.

• \(C^0\): coincide \(f\). Puedes dibujar la curva sin levantar el lápiz. En el knot todavía puede haber una esquina (la pendiente cambia de golpe).
• \(C^1\): coinciden \(f\) y \(f'\). Ya no hay esquina: el flujo y el ritmo al que sube o baja se entregan juntos.
• \(C^2\): coinciden \(f\), \(f'\) y \(f''\). Tampoco hay quiebre de curvatura: la curva no cambia de carácter de golpe al cruzar el corte.

El spline cúbico es \(C^2\). En el tráfico, a las 8:00, la madrugada le entrega a la jornada el mismo caudal, la misma pendiente y el mismo arqueo. En la curva estimada el knot conserva valor, pendiente y segunda derivada.

“Costura” no es el nombre técnico: sirve como imagen. En la literatura se habla de condiciones de continuidad o de empalme en los knots, y de condiciones de frontera en los extremos del intervalo.

Acabas de ver por qué hace falta: partir el día era la idea correcta, pero estimar cada pieza por su lado produce un salto en el knot.

Por qué se llama cúbico

El apellido no viene de que haya tres knots ni tres modelos. Viene del grado de cada pieza.

Recuerda la familia:

• grado 1: una recta, \(p(x)=\beta_0+\beta_1 x\)
• grado 2: una parábola, \(p(x)=\beta_0+\beta_1 x+\beta_2 x^2\)
• grado 3: un cúbico, \(p(x)=\beta_0+\beta_1 x+\beta_2 x^2+\beta_3 x^3\)

Un spline cúbico es esto: en cada región vive un polinomio de grado 3. En el tráfico, con \(\xi_1=8\) y \(\xi_2=17\), escribes primero las tres piezas y, de una vez, sus derivadas. Antes de imponer continuidad o condiciones de frontera, las tres piezas tienen doce coeficientes libres.

Las tres piezas: \(f\), \(f'\) y \(f''\)

\[ \begin{aligned} p_1(x) &= a_1+b_1 x+c_1 x^2+d_1 x^3, && 0\leq x<8,\\ p_2(x) &= a_2+b_2 x+c_2 x^2+d_2 x^3, && 8\leq x<17,\\ p_3(x) &= a_3+b_3 x+c_3 x^2+d_3 x^3, && 17\leq x\leq 24. \end{aligned} \]

La curva del día es \(f(x)=p_i(x)\) en cada tramo. Derivando término a término, la misma receta en las tres regiones, obtienes la pendiente y la curvatura:

\[ p_i'(x)=b_i+2c_i x+3d_i x^2, \qquad p_i''(x)=2c_i+6d_i x, \qquad i=1,2,3. \]

En el tráfico: \(f\) es el flujo, \(f'\) cuánto sube o baja por hora, \(f''\) cómo se arquea esa pendiente. Cada tramo puede curvarse (tiene \(x^2\) y \(x^3\)). Si quedaran independientes, volverías al salto. El spline no las deja independientes: primero simplifica los extremos, después cose los knots.

Primero: condiciones de frontera

En \(x=0\) y \(x=24\) no hay otro tramo con el que empalmar. El spline cúbico natural (una de las variantes más usadas) pide curvatura nula en los extremos:

\[ f''(0)=0, \qquad f''(24)=0. \]

Sustituye en las segundas derivadas que ya tienes:

\[ p_1''(0)=2c_1=0 \implies c_1=0, \]

\[ p_3''(24)=2c_3+6d_3\cdot 24=2c_3+144d_3=0 \implies c_3=-72d_3. \]

Esas dos igualdades reducen las piezas de los extremos. El tramo del medio no toca la frontera, así que sigue con sus cuatro coeficientes. Las piezas de los extremos quedan expresadas mediante los coeficientes restantes:

\[ \begin{aligned} p_1(x) &= a_1+b_1 x+d_1 x^3,\\ p_1'(x) &= b_1+3d_1 x^2,\\ p_1''(x) &= 6d_1 x, \end{aligned} \qquad \begin{aligned} p_3(x) &= a_3+b_3 x-72d_3 x^2+d_3 x^3,\\ p_3'(x) &= b_3-144d_3 x+3d_3 x^2,\\ p_3''(x) &= 6d_3(x-24). \end{aligned} \]

Lee las curvaturas simplificadas: \(p_1''(x)=6d_1 x\) empieza en cero a las 0:00 y cambia de forma lineal con \(x\) (crece en valor absoluto, con el signo de \(d_1\)); \(p_3''(x)=6d_3(x-24)\) llega a cero a las 24:00. No le pides un pico inventado a los extremos; el doblez extra se reserva para los knots del medio, donde el fenómeno sí cambia de régimen.

De doce coeficientes pasas a diez: se fueron \(c_1\) y \(c_3\).

Después: condiciones de continuidad

En cada knot el cúbico de la izquierda y el de la derecha tienen que coincidir en \(f\), \(f'\) y \(f''\). En el tráfico eso se lee en vehículos por minuto, no en abstracto.

Empalme de las 8, la madrugada le entrega el flujo a la jornada:

\[ \begin{aligned} p_1(8) &= p_2(8) && \text{mismo flujo a las 8:00} \\[4pt] p_1'(8) &= p_2'(8) && \text{el flujo sube o baja al mismo ritmo} \\[4pt] p_1''(8) &= p_2''(8) && \text{el flujo se arquea igual: no hay quiebre} \end{aligned} \]

Con \(p_1\) ya simplificado (\(c_1=0\)):

\[ \begin{aligned} a_1+8b_1+512d_1 &= a_2+8b_2+64c_2+512d_2,\\ b_1+192d_1 &= b_2+16c_2+192d_2,\\ 48d_1 &= 2c_2+48d_2. \end{aligned} \]

Si \(p_1(8)\neq p_2(8)\), a las 7:59 predices un caudal y a las 8:01 otro: un salto imposible. Si las pendientes no coinciden, el gráfico forma una esquina. Si fallan las curvaturas, la curva “cambia de carácter” de golpe al cruzar las 8.

Empalme de las 17, la jornada le entrega el flujo a la noche:

\[ \begin{aligned} p_2(17) &= p_3(17) && \text{mismo flujo a las 17:00} \\[4pt] p_2'(17) &= p_3'(17) && \text{misma pendiente} \\[4pt] p_2''(17) &= p_3''(17) && \text{misma curvatura} \end{aligned} \]

Con \(p_3\) ya simplificado (\(c_3=-72d_3\)):

\[ \begin{aligned} a_2+17b_2+289c_2+4913d_2 &= a_3+17b_3-15895d_3,\\ b_2+34c_2+867d_2 &= b_3-1581d_3,\\ 2c_2+102d_2 &= -42d_3. \end{aligned} \]

Seis ecuaciones lineales sobre diez coeficientes. Quedan 4 grados de libertad: esos son los que se ajustan a la nube del tráfico por mínimos cuadrados. El resto del vector no es libre: lo amarran la frontera y la continuidad.

Fíjate en los números del empalme de las 17 (\(15895\), \(1581\)): salen de evaluar potencias \(17^2\), \(17^3\). Esa es la base monomial, y se pone fea rápido. Por eso, en la práctica, no armas esta matriz a mano. SplineTransformer fabrica columnas que ya cumplen las condiciones de continuidad en los knots (no las de frontera del spline natural, como verás enseguida). En la siguiente celda solo eliges knots y grado, y LinearRegression estima los pesos.

Un spline lineal (grado 1 en cada tramo) solo pide continuidad de \(f\): puedes dibujarlo sin levantar el lápiz, pero en el knot puede quebrar. El cúbico pide \(f\), \(f'\) y \(f''\). Por eso se siente de una pieza aunque por debajo sean varios polinomios.

¿Por qué grado 3 y no 10? Un cúbico tiene cuatro números por región. Las condiciones de continuidad y de frontera gastan parte de ellos y dejan libertad para ajustar los datos sin volver a las ondulaciones del polinomio global. El cúbico es el equilibrio que buscas en este capítulo: flexible por regiones, \(C^2\) en los empalmes.

Un poco más abajo vas a ver las tres geometrías juntas: salto, esquina o curva que se desliza. La idea que te tienes que llevar: flexibilidad local con continuidad \(C^2\).

Montar el spline cúbico, paso a paso

Sí: no hay que armar las columnas a mano. scikit-learn ya las calcula con SplineTransformer. Tú eliges el grado y los knots; él fabrica una base local que ya nace \(C^2\). Después, LinearRegression estima los pesos de esas columnas de la base para aproximar el flujo.

knots incluye los extremos del día. En el tráfico: \(0\), \(8\), \(17\) y \(24\). include_bias=False deja el intercepto en el modelo, como en el resto del libro.

Vas a ver 5 columnas. Con el intercepto que agrega LinearRegression son 6 parámetros, no los 4 grados de libertad del spline natural que acabas de contar a mano (esos 4 ya incluían el intercepto: \(a_1\) estaba entre los diez coeficientes). La diferencia son justo las dos condiciones de frontera: SplineTransformer no impone \(f''(0)=f''(24)=0\), así que el conteo es \(12 - 6 = 6\) (doce coeficientes menos seis igualdades de continuidad). Lo que sí heredas es el empalme \(C^2\) en los knots.

bases = SplineTransformer(
    degree=3,  # cúbico: cada pieza es grado 3, empalme C²
    knots=np.array([[0.0], [8.0], [17.0], [24.0]]),  # extremos del día + cortes 8 y 17
    include_bias=False,  # el intercepto lo pone LinearRegression, no la base
)

X_spline = bases.fit_transform(X)  # fabrica las columnas B_j(hora) con los knots que le diste
print(f"Columnas que fabricó SplineTransformer: {X_spline.shape[1]}")  # cuántas bases salieron

# Arriba solo se construyen las bases del spline. 

# abajo ajustamos los pesos

modelo_cubico = LinearRegression()  # la misma recta de siempre, ahora sobre las B_j
modelo_cubico.fit(X_spline, y)  # estima los pesos β_j

# El modelo no come horas crudas: come las columnas B_j. X_grid son horas de la
# grilla (0.57, 0.63, …). transform las convierte en las MISMAS bases, ya evaluadas
# ahí. Si hicieras modelo_cubico.predict(X_grid) fallaría: espera 5 columnas, no 1.
X_spline_grid = bases.transform(X_grid)
pred_cubico_paso_a_paso = modelo_cubico.predict(X_spline_grid)

figura_curvas(
    hora,  # eje X de las observaciones
    y,  # flujo observado
    x_grid,  # grilla donde se dibuja la curva
    curvas={"Spline cúbico": pred_cubico_paso_a_paso},  # predicción ya calculada
    y_verdad=y_real_grid,  # relación verdadera, punteada
    knots=knots_demo,  # verticales en 8 y 17
    colores=[COLOR_SPLINE],
    titulo="SplineTransformer arma las bases; LinearRegression estima los pesos",
)

Columnas que fabricó SplineTransformer: 5

Ahora sí: comparar salto, esquina y curva suave

Ya montaste el cúbico con SplineTransformer. La pregunta es qué ganas frente a lo de antes.

Tres ajustes, mismos datos, mismos cortes \(8\) y \(17\):

• izquierda: las tres rectas sueltas (pueden saltar);
• centro: el mismo SplineTransformer, ahora con degree=1 (continuo, puede tener esquina);
• derecha: el cúbico de la celda anterior (degree=3).

No hay función nueva. Cambia solo el grado.

pred_salto = pred_tramos_grid  # las tres rectas sueltas de antes
pred_cubico_spline = pred_cubico_paso_a_paso  # el cúbico que acabas de ajustar

# Mismos knots, grado 1: continuo, puede quebrar (esquina)
bases_lin = SplineTransformer(
    degree=1,
    knots=np.array([[0.0], [8.0], [17.0], [24.0]]),
    include_bias=False,
)
modelo_lin_sp = LinearRegression()
modelo_lin_sp.fit(bases_lin.fit_transform(X), y)  # ajusta β sobre las bases lineales
pred_lineal_spline = modelo_lin_sp.predict(bases_lin.transform(X_grid))  # predice en la grilla

# ¿Cuánto cuesta la continuidad con solo dos knots? RMSE in-sample de los tres ajustes
rmse_tramos = np.sqrt(mean_squared_error(y, pred_tramos_obs))
rmse_lin_sp = np.sqrt(mean_squared_error(y, modelo_lin_sp.predict(bases_lin.transform(X))))
rmse_cub_sp = np.sqrt(mean_squared_error(y, modelo_cubico.predict(bases.transform(X))))
print(f"RMSE in-sample  |  tramos: {rmse_tramos:.3f}  |  spline lineal: {rmse_lin_sp:.3f}  |  spline cúbico: {rmse_cub_sp:.3f}")
print(f"Mínimo del spline cúbico en la malla: {pred_cubico_spline.min():.1f} vehículos/min")

figura_paneles(
    hora, y, x_grid,
    {
        "Tramos independientes\n(pueden saltar)": pred_salto,
        "Spline lineal\n(continuo, con esquinas)": pred_lineal_spline,
        "Spline cúbico\n(continuo en nivel, pendiente y curvatura)": pred_cubico_spline,
    },
    knots=knots_demo,
    trozos=piezas_lineales,  # para que el salto de la izquierda no se tape
    colores=[COLOR_POLI, COLOR_LINEAL, COLOR_SPLINE],
    titulo_figura="La diferencia no es partir el dominio: es cómo se unen las piezas",
)

RMSE in-sample  |  tramos: 17.556  |  spline lineal: 27.569  |  spline cúbico: 24.761
Mínimo del spline cúbico en la malla: -38.9 vehículos/min

Las verticales punteadas de la figura señalan los cortes de las 8:00 y las 17:00. Entre el panel izquierdo y el central cambia si hay condiciones de continuidad; entre el central y el derecho cambia, además, el grado de las piezas.

A la izquierda, cada recta vive su vida: busca el salto. Al centro, el spline lineal es continuo, lo dibujas sin levantar el lápiz, pero en el knot puede quebrar: solo se impuso \(f\). A la derecha, el cúbico es \(C^2\): mismo flujo, misma pendiente, misma curvatura. Las 5 columnas que fabricó SplineTransformer ya traían esas igualdades; LinearRegression solo estimó los pesos, más el intercepto.

Ojo con una lectura apresurada: más suave no significa mejor ajuste. Con solo dos knots internos, las condiciones de continuidad amarran las piezas y les quitan libertad. El RMSE in-sample lo muestra: 17.556 con las tres rectas sueltas, 27.569 con el spline lineal y 24.761 con el cúbico, que además se hunde hasta −38.9 vehículos por minuto a las 23:25. La continuidad arregla el empalme, no la falta de flexibilidad; para devolverle flexibilidad al spline hay que darle más knots, y eso es lo que viene.

20.9 Qué controla la flexibilidad: el número de knots

En el capítulo de regularización el parámetro que controlaba la complejidad era \(\alpha\). Aquí el control más directo es cuántos knots permites.

Pocos knots  →  modelo rígido     →  alto sesgo (underfitting)
Muchos knots →  modelo flexible   →  alta varianza (overfitting)

El mecanismo es otro; el problema de fondo es el mismo: el compromiso sesgo-varianza (bias-variance tradeoff), es decir, controlar la complejidad.

Ahora ajustas splines cúbicos con 3, 5, 8 y 15 knots internos, colocados de forma uniforme. No elijas todavía el “mejor”: mira qué le ocurre a la curva.

El mismo spline cúbico, cuatro cantidades de knots

Cuatro paneles, el mismo dataset. 3, 5, 8 y 15 knots internos. Los cortes se marcan para que veas dónde se le permitió cambiar al modelo.

Al algoritmo no le pasas horas elegidas a simple vista. Le pasas cuántos knots y una regla. Con knots="uniform" (el valor por defecto), fit(X) reparte ese número entre \(x_{\min}\) y \(x_{\max}\) en intervalos iguales: es el mismo np.linspace de la celda de abajo, que vuelve a aparecer en la sección sobre dónde colocar los knots. En esa celda lo calculas a mano y se lo entregas en knots=..., para que veas las posiciones. Las verticales de la figura son esas mismas horas, no otro cálculo paralelo.

n_knots_internos = [3, 5, 8, 15]
paneles = {}
cortes = []
filas_flex = []  # una fila por cantidad de knots: MAE, RMSE, R² in-sample
xmin, xmax = float(hora.min()), float(hora.max())  # el rango de X: ahí se reparte el reloj

for n_k in n_knots_internos:
    # K internos + 2 extremos, a la misma distancia en horas. Esto SÍ se lo pasamos
    # al transformador: no adivina las 8 y las 17; las calcula así, sobre [xmin, xmax].
    knots_todos = np.linspace(xmin, xmax, n_k + 2).reshape(-1, 1)
    knots_internos = knots_todos[1:-1, 0]  # las verticales de la figura: sin xmin ni xmax

    bases_k = SplineTransformer(
        degree=3,
        knots=knots_todos,  # posiciones, no solo el número; incluye los extremos
        include_bias=False,
    )
    # Equivalente a n_knots=n_k+2 y knots="uniform": fit(X) haría este linspace solo.
    modelo_k = LinearRegression()
    modelo_k.fit(bases_k.fit_transform(X), y)
    pred = modelo_k.predict(bases_k.transform(X_grid))
    pred_obs = modelo_k.predict(bases_k.transform(X))
    mae, rmse, r2 = metricas_regresion(y, pred_obs, f"Spline cúbico, {n_k} knots internos (in-sample)")
    paneles[f"{n_k} knots internos  |  RMSE={rmse:.1f}  $R^2$={r2:.2f}"] = pred
    cortes.append(knots_internos)  # las mismas horas que recibió SplineTransformer
    filas_flex.append(
        {
            "Knots internos": n_k,
            "MAE": mae,
            "RMSE": rmse,
            "R²": r2,
        }
    )

figura_paneles(
    hora, y, x_grid, paneles,
    y_verdad=y_real_grid,
    knots_por_panel=cortes,
    colores=[COLOR_SPLINE] * 4,
    n_cols=2,
    titulo_figura="Pocos knots: rígido. Más knots: la curva alcanza los dos picos",
)

# Mismas métricas de siempre, ahora en una tabla. In-sample: el error es sobre
# los puntos con los que se ajustó, no sobre un día nuevo.
tabla_flex = pd.DataFrame(filas_flex)
print("Métricas in-sample (mismo martes con el que se ajustó)")
tabla_flex.round(3)

	Knots internos	MAE	RMSE	R²
0	3	15.559	18.923	0.632
1	5	10.993	13.510	0.812
2	8	5.972	7.501	0.942
3	15	5.036	6.264	0.960

Spline cúbico, 3 knots internos (in-sample)
  MAE  = 15.559
  RMSE = 18.923
  R²   = 0.632
Spline cúbico, 5 knots internos (in-sample)
  MAE  = 10.993
  RMSE = 13.510
  R²   = 0.812
Spline cúbico, 8 knots internos (in-sample)
  MAE  = 5.972
  RMSE = 7.501
  R²   = 0.942
Spline cúbico, 15 knots internos (in-sample)
  MAE  = 5.036
  RMSE = 6.264
  R²   = 0.960
Métricas in-sample (mismo martes con el que se ajustó)

En la tabla, de 3 a 15 knots internos, MAE y RMSE bajan y \(R^2\) sube: 15 knots gana las tres columnas (MAE 5.036, RMSE 6.264, \(R^2\) 0.960). A más cortes, más permiso para pegarse a estos puntos. Eso no elige modelo: las métricas son in-sample, medidas sobre el mismo martes con el que se ajustó, ruido incluido. La selección requiere estimar el error sobre observaciones reservadas del mismo fenómeno.

En la figura, con 3 knots internos el modelo apenas puede cambiar de pieza: deja los picos en ~78 y ~81 cuando la curva verdadera llega a ~109 y ~114, y hacia las 2:30 se hunde hasta ~3. Eso es subajuste (underfitting, alto sesgo). Con 5 y 8 la curva empieza a parecerse al día real. Con 15 knots internos sigue de cerca los dos picos. Como aquí conoces la función verdadera, puedes medir su distancia a ella: 2.7 (RMSE sobre la malla) frente a 5.3 con 8 knots y 18.3 con 3. En esta muestra es la mejor aproximación de las cuatro, lo que no excluye que una parte de su ajuste sea ruido. Ese es el punto delicado: en la tabla in-sample, 15 knots ganaría aunque estuviera sobreajustando, porque el error de entrenamiento casi siempre baja al agregar knots. Con muchos más knots aparecería el sobreajuste (overfitting, alta varianza): ondulaciones que son ruido de este martes, aunque la función siga viéndose “suave” comparada con un polinomio de grado 12. Para encontrar la frontera necesitas la validación cruzada de más abajo; la tabla de entrenamiento no la muestra.

El espectro es el de siempre:

subajuste \(\leftrightarrow\) complejidad adecuada \(\leftrightarrow\) sobreajuste

Ridge y Lasso penalizan el tamaño de los coeficientes (\(\alpha\) hace el papel del \(\lambda\) de otros textos). Aquí el número de knots controla la flexibilidad local. Distinto mecanismo, mismo compromiso entre sesgo y varianza. El equivalente de «meter todas las variables y no penalizar» sería aquí poner muchísimos knots sin validar. Y, igual que entonces, el ajuste sobre los datos usados para entrenar no es el juez.

20.10 Dónde colocar los knots

Dos estrategias simples bastan para este capítulo. Los knots uniformes en \(X\) reparten el rango observado en intervalos de igual longitud; es una regla transparente que funciona cuando las observaciones cubren el eje de forma más o menos pareja. Los knots por cuantiles de \(X\) se colocan en percentiles de la muestra (terciles, cuartiles); si \(X\) está amontonada en una región, con muchas mediciones en hora pico y casi ninguna a las 4 a.m., los cuantiles gastan más knots donde hay más datos, que es donde puedes permitirte estimar cambios locales.

Ninguna de las dos reglas es preferible en todos los casos. El error que sí debes evitar es de procedimiento: fijar los knots inspeccionando \(y\), por ejemplo colocándolos en los picos del tráfico, sobre el conjunto completo, y evaluar luego el modelo como si esas posiciones fueran independientes de la respuesta. En ese caso \(y\) interviene dos veces, primero en el diseño del spline y después en su evaluación. La salida es usar una regla que no mire \(y\) y calcularla solo con cada partición de entrenamiento de la validación cruzada; el Pipeline se encarga de esto último, pero no te protege si escribes a mano knots elegidos mirando \(y\).

Un eje X desigual: uniformes frente a cuantiles

En este ejemplo no ajustas ningún modelo: no hay variable respuesta ni mínimos cuadrados, y solo ves dónde caen los knots cuando las horas no están repartidas de forma pareja. Fabricas un día artificial con muchas mediciones en las dos horas pico y pocas de madrugada y de noche, y sobre esas horas colocas \(K=6\) knots internos con las dos reglas: los uniformes quedan a la misma distancia en el reloj aunque en un tramo haya 55 puntos y en otro 4, y los de cuantiles dejan (casi) el mismo número de observaciones entre knot y knot.

Si las marcas coinciden, el muestreo era parejo y da lo mismo qué regla uses. Si se separan, sobre todo cerca de las 8 y de las 18, estás viendo el efecto de la densidad de \(X\) sobre el diseño del spline.

# Un día mal muestreado a propósito: muchas mediciones en las dos horas pico,
# pocas de madrugada y de noche. Así se ve la diferencia entre las dos reglas.
hora_desigual = np.sort(  # de menor a mayor: el histograma se lee como un eje de tiempo
    np.concatenate(  # pega los cinco bloques en un solo vector de horas
        [
            rng.uniform(0.5, 6.5, 20),    # madrugada: 20 puntos en ~6 h
            rng.uniform(7.0, 9.5, 55),    # pico mañana: 55 puntos en ~2.5 h
            rng.uniform(10.0, 16.0, 25),  # mediodía: 25 puntos en 6 h
            rng.uniform(16.5, 19.5, 55),  # pico tarde: 55 puntos en ~3 h
            rng.uniform(20.0, 23.5, 20),  # noche: 20 puntos en ~3.5 h
        ]
    )
)

n_internos = 6  # K knots INTERIORES; los extremos no se cuentan (el spline ya los usa)

# Uniformes: K+2 puntos igualmente espaciados entre xmin y xmax (incluidos).
# [1:-1] tira xmin y xmax. El paso es (xmax-xmin)/(K+1): reparte el RELOJ.
knots_unif = np.linspace(hora_desigual.min(), hora_desigual.max(), n_internos + 2)[1:-1]

# Cuantiles: las mismas K+2 probabilidades de 0 a 1; se tiran 0 y 1 (min y max
# muestrales) y se evalúa el percentil de X. Equidistancia en el ORDEN de la
# muestra: entre knot y knot hay (casi) el mismo número de observaciones, no las mismas horas.
knots_cuant = np.quantile(hora_desigual, np.linspace(0, 1, n_internos + 2)[1:-1])

figura_histograma_knots(
    hora_desigual,  # histograma de X arriba; marcas de knots abajo
    {
        "Uniformes": (knots_unif, COLOR_LINEAL, "--"),  # espacio constante en horas
        "Cuantiles": (knots_cuant, COLOR_SPLINE, "-"),  # espacio constante en cuantiles
    },
    titulo_hist="X no es uniforme: hay muchas más observaciones en las horas pico",
    titulo_knots="Los cuantiles agrupan knots donde hay datos; los uniformes, donde hay rango",
)

print("Knots uniformes :", np.round(knots_unif, 2))  # deberían verse casi aritméticos
print("Knots cuantiles :", np.round(knots_cuant, 2))  # más apretados cerca de 8 y 18

Knots uniformes : [ 4.39  7.57 10.76 13.94 17.12 20.3 ]
Knots cuantiles : [ 7.13  8.61  9.93 16.17 17.74 18.99]

El histograma de arriba no es plano: la ciudad midió más en las horas críticas. Los cortes rojos (uniformes) se espacian por reloj: de 4.39 a 20.30 h, cada 3.2 h. Los verdes (cuantiles) se concentran en los dos bloques pico: tres entre 7.13 y 9.93 h y tres entre 16.17 y 18.99 h. No ponen ninguno en la madrugada, el mediodía ni la noche. Si tu interés es estimar bien los picos y allí es donde tienes muestra, los cuantiles son una regla razonable. Si tu interés es cubrir el reloj de forma equitativa aunque a las 4 a.m. tengas cuatro puntos, los uniformes cubren mejor el dominio… y son más frágiles en esas zonas vacías.

En la siguiente sección esas dos reglas son las opciones knots="uniform" (la que usa el capítulo, y la de por defecto) y knots="quantile" del transformador. Aquí basta la geometría: dónde hay datos, ahí puedes permitirte un knot.

20.11 Implementación con scikit-learn: SplineTransformer

Ya lo usaste a mano: knots, grado y LinearRegression. Ahora lo dejas en un Pipeline. Es la forma que vas a usar cuando elijas n_knots con k-fold, y cuando las dos reglas de colocación (uniforme o cuantiles) tengan que calcularse dentro de cada fold, no con el dataset entero.

El spline no reemplaza la regresión lineal. Transforma \(X\) en nuevas variables y después una regresión lineal aprende los pesos. Esta construcción es una transformación de características: convierte la hora en columnas sobre las que se ajusta una regresión.

X
 ↓
SplineTransformer
 ↓
B₁(X), B₂(X), ..., Bₘ(X)
 ↓
LinearRegression
 ↓
ŷ

n_knots en scikit-learn cuenta también los knots de los extremos. degree=3 es el spline cúbico. include_bias=False deja el intercepto en manos de LinearRegression. Las dos reglas de la sección anterior son knots="uniform" (por defecto) y knots="quantile".

Qué produjo realmente SplineTransformer

Antes del Pipeline, mira la transformación sola: dimensión de \(X\) frente a dimensión de las bases, y unas pocas filas. Esas columnas son las \(B_j(x)\). La regresión posterior ajusta sus coeficientes sobre las columnas transformadas, manteniendo la partición de entrenamiento.

spline_feat = SplineTransformer(n_knots=5, degree=3, include_bias=False)
X_spline = spline_feat.fit_transform(X)

print(f"X.shape        = {X.shape}")
print(f"X_spline.shape = {X_spline.shape}")
print("Primeras filas de las bases B_j(X):")
print(
    pd.DataFrame(
        X_spline[:6],
        columns=[f"B{j+1}" for j in range(X_spline.shape[1])],
    ).round(3).to_string(index=False)
)

pipeline_spline = Pipeline(
    [
        ("spline", SplineTransformer(n_knots=5, degree=3, include_bias=False)),
        ("modelo", LinearRegression()),
    ]
)
pipeline_spline.fit(X, y)
pred_sk_obs = pipeline_spline.predict(X)
pred_sk_grid = pipeline_spline.predict(X_grid)

print(
    f"Columnas del spline en el Pipeline: {pipeline_spline.named_steps['modelo'].coef_.shape[0]}"
)
_ = metricas_regresion(y, pred_sk_obs, "Pipeline SplineTransformer(n_knots=5, degree=3)")

figura_curvas(
    hora, y, x_grid,
    curvas={"Spline (Pipeline)": pred_sk_grid},
    y_verdad=y_real_grid,
    colores=[COLOR_SPLINE],
    titulo="SplineTransformer expande X; LinearRegression estima los pesos",
)

X.shape        = (180, 1)
X_spline.shape = (180, 6)
Primeras filas de las bases B_j(X):
   B1    B2    B3  B4  B5  B6
0.167 0.667 0.167 0.0 0.0 0.0
0.139 0.663 0.197 0.0 0.0 0.0
0.137 0.663 0.200 0.0 0.0 0.0
0.137 0.663 0.200 0.0 0.0 0.0
0.123 0.658 0.218 0.0 0.0 0.0
0.113 0.653 0.234 0.0 0.0 0.0
Columnas del spline en el Pipeline: 6
Pipeline SplineTransformer(n_knots=5, degree=3)
  MAE  = 15.559
  RMSE = 18.923
  R²   = 0.632

\(X\) tenía una columna; X_spline tiene seis: n_knots + degree − 1 = 7 bases, menos la que quita include_bias=False. Cada una es una base \(B_j\) evaluada en cada hora. Eso es la ingeniería de variables (feature engineering): no cambiaste el modelo, cambiaste las variables.

Fíjate en las métricas: RMSE 18.923 y \(R^2\) 0.632, idénticas al panel de “3 knots internos” de la figura de flexibilidad. No es casualidad: n_knots=5 cuenta los dos extremos, así que es el mismo modelo, y por eso la curva todavía se queda corta en los dos picos. Lo que cambió es la forma de escribirlo, no la flexibilidad.

El Pipeline hace exactamente esa expansión y después LinearRegression estima los \(\beta_j\). Si más adelante cambias n_knots o degree, no cambias de algoritmo: cambias el espacio de variables. Por eso el spline se lleva tan bien con pipelines, k-fold y métricas. Y por eso la transformación tiene que vivir dentro del Pipeline: si no, los folds de prueba ya vieron cómo se construyeron las columnas.

20.12 Elegir la flexibilidad con datos independientes de las ilustraciones

Las figuras anteriores explican las formas posibles con una primera muestra de 180 mediciones. Para seleccionar y evaluar usas la segunda muestra independiente, cuya prueba se reservó en la celda de generación. Hay 135 observaciones de entrenamiento y 45 de prueba. Todos los candidatos ajustan transformaciones y coeficientes dentro de los folds de entrenamiento.

Los repartos aleatorios miden nuevas observaciones del mismo mecanismo sintético. Pronosticar otro día exigiría mediciones de varios días y una validación temporal. Las alusiones al tráfico ayudan a interpretar unidades, pero no convierten esta simulación en evidencia de pronóstico diario.

Recorrer n_knots con k-fold

En esta celda eliges cuántos knots usar con validación cruzada, para no quedarte con el modelo que mejor se ve en entrenamiento. Pruebas una lista de candidatos, de pocos knots a 40, para ver el descenso del RMSE y, si la flexibilidad sobra, su repunte a la derecha; más allá de 40, cada fold de entrenamiento tiene apenas unos 108 puntos, algunos tramos entre knots quedan casi vacíos y el ajuste se vuelve numéricamente inestable. Para cada valor armas el mismo Pipeline (SplineTransformer y LinearRegression) y estimas el RMSE solo con el entrenamiento, con el mismo KFold de 5 folds, y conservas el n_knots con menor RMSE promedio.

En scikit-learn, n_knots cuenta también los extremos \(x_{\min}\) y \(x_{\max}\), así que los knots internos son \(n_{\mathrm{knots}}-2\). Por ejemplo, n_knots=5 son 3 knots internos, no el panel de “5 knots internos” de la figura de flexibilidad, que equivaldría a n_knots=7; y el panel de “15 knots internos” equivale a n_knots=17.

# La prueba de la muestra independiente ya se reservó al generar los datos.
print(f"Entrenamiento: {len(y_ent)} observaciones  |  prueba: {len(y_prueba)}")

kf = KFold(n_splits=5, shuffle=True, random_state=SEMILLA)
valores_knots = [3, 4, 5, 6, 8, 10, 12, 15, 20, 25, 30, 40]

filas_knots = []
for n_knots in valores_knots:
    pipe = Pipeline(
        [
            ("spline", SplineTransformer(n_knots=n_knots, degree=3, include_bias=False)),
            ("modelo", LinearRegression()),
        ]
    )
    rmse_folds = -cross_val_score(
        pipe, X_ent, y_ent, cv=kf, scoring="neg_root_mean_squared_error"
    )
    filas_knots.append(
        {
            "n_knots": n_knots,
            "RMSE_CV": rmse_folds.mean(),
            "STD_RMSE": rmse_folds.std(ddof=1),
        }
    )

tabla_knots = pd.DataFrame(filas_knots)
mejor_n_knots = int(tabla_knots.loc[tabla_knots["RMSE_CV"].idxmin(), "n_knots"])
print(f"Mejor n_knots según RMSE de CV: {mejor_n_knots}")
tabla_knots.round(3)

	n_knots	RMSE_CV	STD_RMSE
0	3	25.539	2.216
1	4	24.965	2.188
2	5	20.019	1.740
3	6	15.941	2.231
4	8	14.829	0.956
5	10	8.228	0.472
6	12	9.665	0.519
7	15	7.477	1.208
8	20	6.801	0.665
9	25	6.594	0.516
10	30	6.595	0.434
11	40	7.713	0.470

Entrenamiento: 135 observaciones  |  prueba: 45
Mejor n_knots según RMSE de CV: 25

Complejidad frente a error de validación

Cada punto es el RMSE promedio de los 5 folds. Las barras indican la desviación estándar de ese RMSE entre folds: miden cuánto varía el error al cambiar la partición, no un intervalo de confianza. Si dos valores de n_knots tienen medias parecidas y las barras se solapan, tómalo como una señal de incertidumbre, no como una prueba: las barras describen la dispersión entre folds, y todos los candidatos se evalúan sobre los mismos folds.

figura_cv(
    tabla_knots["n_knots"].to_list(),
    tabla_knots["RMSE_CV"].to_list(),
    tabla_knots["STD_RMSE"].to_list(),
    elegido=mejor_n_knots,
    titulo="La complejidad ahora se elige con k-fold, no con el gráfico de entrenamiento",
)

En la muestra independiente de selección, el mínimo de la rejilla está en n_knots=25, con RMSE CV de 6.594 y desviación entre folds de 0.516. Estos knots incluyen los dos extremos; no son el número de cortes interiores de las figuras ilustrativas.

La regla heurística de un error estándar usaría 0.516/√5 ≈ 0.231 y un umbral de 6.825; entre los candidatos que lo cumplen, el menor n_knots es 20. Esa regla favorece sencillez y no es una prueba estadística: los folds comparten entrenamiento. Aquí conservas el mínimo predefinido para el reporte.

Los RMSE de la rejilla no tienen que ser monótonos: cambian la ubicación de las bases y la sensibilidad a cada fold. El mínimo seleccionado puede ser optimista, por lo que falta abrir la prueba que quedó reservada.

20.13 Comparación final: lineal, polinómico y spline

En esta sección comparas tres familias (recta, polinomio y spline cúbico) para decidir cuál generaliza mejor en el tráfico, con el mismo KFold de 5 folds y las mismas métricas: MAE, RMSE y \(R^2\) de validación cruzada, con su promedio y su desviación entre folds.

La recta no tiene hiperparámetros y entra tal cual. Al polinomio le eliges el grado por RMSE de validación cruzada sobre el entrenamiento, y el spline ya trae su n_knots de la celda anterior. Los tres se evalúan después en la misma tabla, primero con validación cruzada sobre el entrenamiento y, al final, una sola vez sobre la prueba. El mismo protocolo para los tres evita que la comparación dependa de haber favorecido a uno a simple vista, y pone a prueba una idea: el modelo más flexible a la vista no es necesariamente el que mejor generaliza.

Elegir el grado del polinomio con el mismo k-fold

Grados 2, 3, 5, 8, 10, 12, 14, 16, 18 y 20. El ganador entra a la tabla final; los demás no se esconden: ves su RMSE CV.

candidatos_grado = [2, 3, 5, 8, 10, 12, 14, 16, 18, 20]
filas_grado = []
for grado in candidatos_grado:
    pipe = pipeline_polinomico(grado)
    rmse_folds = -cross_val_score(
        pipe, X_ent, y_ent, cv=kf, scoring="neg_root_mean_squared_error"
    )
    filas_grado.append(
        {
            "grado": grado,
            "RMSE_CV_promedio": rmse_folds.mean(),
            "RMSE_CV_std": rmse_folds.std(ddof=1),
        }
    )

tabla_grado = pd.DataFrame(filas_grado)
mejor_grado = int(tabla_grado.loc[tabla_grado["RMSE_CV_promedio"].idxmin(), "grado"])
print(f"Mejor grado polinómico según RMSE de CV: {mejor_grado}")
tabla_grado.round(3)

assert candidatos_grado[0] < mejor_grado < candidatos_grado[-1], "Amplía la rejilla de grados antes de evaluar la prueba."

Mejor grado polinómico según RMSE de CV: 14

Tabla de generalización de los tres modelos

Tres pipelines y las tres métricas del libro (MAE, RMSE y \(R^2\)), promediadas entre los cinco folds del entrenamiento; para el RMSE se muestra además la desviación entre folds. La última columna es el RMSE sobre el conjunto de prueba.

def evaluar_cv(modelo, X, y, kf):
    mae = -cross_val_score(modelo, X, y, cv=kf, scoring="neg_mean_absolute_error")
    rmse = -cross_val_score(modelo, X, y, cv=kf, scoring="neg_root_mean_squared_error")
    r2 = cross_val_score(modelo, X, y, cv=kf, scoring="r2")
    return {
        "MAE_CV": mae.mean(),
        "MAE_std": mae.std(ddof=1),
        "RMSE_CV": rmse.mean(),
        "RMSE_std": rmse.std(ddof=1),
        "R2_CV": r2.mean(),
        "R2_std": r2.std(ddof=1),
    }


modelos_finales = {
    "Regresión lineal": LinearRegression(),
    f"Polinomio grado {mejor_grado}": pipeline_polinomico(mejor_grado),
    f"Spline cúbico (n_knots={mejor_n_knots})": Pipeline(
        [
            (
                "spline",
                SplineTransformer(n_knots=mejor_n_knots, degree=3, include_bias=False),
            ),
            ("modelo", LinearRegression()),
        ]
    ),
}

filas_final = []
for nombre, modelo in modelos_finales.items():
    res = evaluar_cv(modelo, X_ent, y_ent, kf)  # CV solo sobre el entrenamiento
    modelo.fit(X_ent, y_ent)  # ajuste final con todo el entrenamiento
    pred_prueba = modelo.predict(X_prueba)  # la prueba se usa una sola vez, aquí
    res["RMSE_prueba"] = float(np.sqrt(mean_squared_error(y_prueba, pred_prueba)))
    filas_final.append({"Modelo": nombre, **res})

tabla_final = pd.DataFrame(filas_final)
tabla_mostrar = pd.DataFrame(
    {
        "Modelo": tabla_final["Modelo"],
        "MAE CV": tabla_final["MAE_CV"].round(3),
        "RMSE CV": [
            f"{m:.3f} ± {s:.3f}"
            for m, s in zip(tabla_final["RMSE_CV"], tabla_final["RMSE_std"])
        ],
        "R² CV": tabla_final["R2_CV"].round(3),
        "RMSE prueba": tabla_final["RMSE_prueba"].round(3),
    }
)
tabla_mostrar

	Modelo	MAE CV	RMSE CV	R² CV	RMSE prueba
0	Regresión lineal	24.189	30.655 ± 2.337	-0.042	24.897
1	Polinomio grado 14	7.015	8.562 ± 0.784	0.917	10.485
2	Spline cúbico (n_knots=25)	5.245	6.594 ± 0.516	0.952	6.306

La recta obtiene RMSE CV de 30.655 y \(R^2\) CV de -0.042. El polinomio elegido tiene grado 14 y RMSE CV de 8.562. El spline con n_knots=25 obtiene 6.594. Las configuraciones polinómica y spline se eligieron con estos folds; sus mínimos describen la selección y pueden ser optimistas.

Los RMSE sobre las 45 observaciones de prueba, en vehículos por minuto, son 24.897 para la recta, 10.485 para el polinomio y 6.306 para el spline. Esta muestra se separó antes de inspeccionarla y es independiente de las ilustraciones iniciales. El reporte confirma el resultado de estos procedimientos ya fijados, sin volver a elegir a partir de la prueba.

El ruido generador tiene desviación 6.5. Es una referencia del error irreducible en esperanza bajo pérdida cuadrática, no una cota de cada RMSE observado. Un resultado de prueba por debajo de 6.5 puede aparecer por azar de las observaciones y no contradice esa referencia.

20.14 Visualización de los modelos evaluados

Las curvas siguientes se ajustan únicamente con las 135 observaciones de entrenamiento de la muestra de evaluación y conservan la configuración elegida por CV. La prueba no se incorpora al dibujo ni al reajuste. La relación generadora se muestra porque es conocida en esta simulación.

pred_lin_grid = LinearRegression().fit(X_ent, y_ent).predict(X_grid)
pred_poli_grid = pipeline_polinomico(mejor_grado).fit(X_ent, y_ent).predict(X_grid)
modelo_spline_final = Pipeline(
    [
        ("spline", SplineTransformer(n_knots=mejor_n_knots, degree=3, include_bias=False)),
        ("modelo", LinearRegression()),
    ]
)
pred_spl_grid = modelo_spline_final.fit(X_ent, y_ent).predict(X_grid)

figura_curvas(
    X_ent.ravel(), y_ent, x_grid,
    {
        "Regresión lineal": pred_lin_grid,
        f"Polinomio grado {mejor_grado}": pred_poli_grid,
        f"Spline cúbico (n_knots={mejor_n_knots})": pred_spl_grid,
    },
    y_verdad=y_real_grid,
    colores=[COLOR_LINEAL, COLOR_POLI, COLOR_SPLINE],
    titulo="Qué gana el spline: picos locales sin comprar ondulaciones globales",
    loc_leyenda="upper left",
)

Las curvas de recta, polinomio y spline se comparan ahora con la misma muestra de entrenamiento del reporte. La recta tiene una sola pendiente; el polinomio usa potencias globales y el spline una base por tramos. Observa dónde cada uno aproxima los picos de las 8:00 y las 18:12 y contrasta esa lectura con los RMSE de la tabla, sin cambiar la configuración después de ver la prueba.

20.15 Errores comunes al trabajar con splines

1. Pensar que más knots siempre es mejor. Más knots es más flexibilidad. La flexibilidad de más memoriza el ruido. El juez es la validación, no el RMSE de entrenamiento ni el “se ve más pegado a los puntos”.

2. Colocar los knots mirando todo el dataset, a veces incluso mirando \(y\), y evaluar después como si no se hubieran usado. Eso es fuga. Los knots (uniformes o por cuantiles) deben nacer del fold de entrenamiento. El Pipeline conserva la misma transformación y el mismo modelo al pasar del entrenamiento a la predicción.

3. Confundir spline con regresión polinómica. Un polinomio de grado 10 es una función global. Un spline cúbico con varios knots es una función por tramos con costura. Pueden tener un número parecido de parámetros y comportamientos opuestos en los extremos.

4. Olvidar la transformación dentro del Pipeline. Ajustar el SplineTransformer a \(X\) completo, transformar, y luego meter las columnas ya listas en cross_val_score es exactamente el error de preprocesar fuera de la validación.

5. Interpretar cada \(\beta_j\) como si fuera una pendiente simple. “El tercer coeficiente es 12, luego cada hora extra suma 12 vehículos” es falso. Cada \(\beta_j\) pondera una columna de SplineTransformer, que actúa sobre un tramo de \(X\), no sobre todo el día. La pendiente del spline cambia a lo largo de \(X\); eso era el punto.

6. Extrapolar demasiado lejos del rango observado. Fuera del intervalo de \(X\) con el que se construyeron las bases, el spline ya no está anclado. No es una ley física: preguntarle el flujo a una hora que no existió en la muestra es pedirle que invente. En SplineTransformer el valor por defecto es extrapolation="constant": fuera del rango de entrenamiento la predicción se queda plana en el valor del borde; con "linear" o "continue" prolonga la tendencia del borde. Ninguna opción convierte la extrapolación en un dato.

7. Evaluar únicamente el ajuste visual. Una curva que recorre los puntos puede ser un sobreajuste elegante. La gráfica convence; el k-fold decide. Usa las dos, en ese orden de autoridad cuando entren en conflicto.

20.16 Dónde aparece esto en inteligencia artificial y analítica de datos

En este capítulo usaste la base spline como una transformación de variables: convierte una columna en varias para que una regresión lineal pueda aprender una curva. Esa idea, la ingeniería de variables, es trabajo diario en analítica y en aprendizaje automático, y los splines aparecen en ella de formas concretas:

• Variables de tiempo en modelos de pronóstico. La hora del día, el día de la semana o el mes son circulares: después de las 23 vienen las 0. Un spline periódico se lo enseña al modelo. Es una práctica habitual al preparar variables para pronosticar demanda, consumo o tráfico, y es el ejemplo que sigue.
• Modelos aditivos interpretables. Un modelo aditivo generalizado (GAM) suma un efecto suave por cada variable, que suele representarse con splines, y lo conecta con la respuesta mediante una función de enlace. Como el efecto de cada variable se puede dibujar por separado, se usa donde el modelo debe poder explicarse, por ejemplo en riesgo de crédito o en salud. Es el tema del capítulo de modelos aditivos.
• Knots elegidos por el algoritmo. El método MARS (multivariate adaptive regression splines) responde de forma automática la pregunta de dónde colocar los knots: prueba cortes, conserva los que reducen el error y después poda los que no compensan su complejidad.
• Redes neuronales con activación ReLU. Una red con ReLU calcula una función lineal por tramos. Con una sola variable de entrada es un spline lineal como el de las esquinas de este capítulo, con knots que no se fijan de antemano: los aprende el entrenamiento. Con varias variables la idea se generaliza a regiones del espacio, cada una con su propia función lineal.

En los cuatro casos la decisión de fondo es la de este capítulo: cuántas piezas permitir, y validarlo con datos que el modelo no vio.

Ejemplo: la hora es circular

El spline que elegiste tiene un defecto que todavía no miraste: trata las 0:00 y las 24:00 como extremos opuestos, cuando son el mismo instante. Medianoche es un empalme que nadie cosió.

SplineTransformer puede coserlo con extrapolation="periodic": obliga a la curva a llegar a las 24:00 con el mismo valor, la misma pendiente y la misma curvatura con que sale de las 0:00. Es conocimiento del problema convertido en una restricción del modelo. Comparas los dos splines con el mismo número de knots, el mismo entrenamiento y el mismo k-fold. No es una comparación que aísle solo la periodicidad: el spline periódico reparte sus knots sobre el día completo, de 0 a 24, y el otro sobre las horas observadas, y por eso terminan con distinto número de columnas.

# La hora es circular: las 23:59 y las 00:01 son vecinas. El spline elegido no lo sabe.
spline_normal = Pipeline(
    [
        ("spline", SplineTransformer(n_knots=mejor_n_knots, degree=3, include_bias=False)),
        ("modelo", LinearRegression()),
    ]
)
spline_periodico = Pipeline(
    [
        (
            "spline",
            SplineTransformer(
                degree=3,
                knots=np.linspace(0, 24, mejor_n_knots).reshape(-1, 1),  # knots sobre el día completo
                extrapolation="periodic",  # lo que sale por las 24:00 entra por las 0:00
                include_bias=False,
            ),
        ),
        ("modelo", LinearRegression()),
    ]
)

filas_periodico = []
for nombre, modelo in {"Spline cúbico": spline_normal, "Spline cúbico periódico": spline_periodico}.items():
    rmse_cv = -cross_val_score(modelo, X_ent, y_ent, cv=kf, scoring="neg_root_mean_squared_error")
    modelo.fit(X_ent, y_ent)
    flujo_24, flujo_0 = modelo.predict([[24.0], [0.0]])  # las 24:00 y las 0:00 son el mismo instante
    filas_periodico.append(
        {
            "Modelo": nombre,
            "Columnas": modelo.named_steps["modelo"].coef_.shape[0],
            "RMSE CV": rmse_cv.mean(),
            "Flujo a las 24:00": flujo_24,
            "Flujo a las 0:00": flujo_0,
            "Diferencia": abs(flujo_24 - flujo_0),
        }
    )

# De las 20:00 a las 04:00 del día siguiente: el eje sigue de largo, el reloj vuelve a 0.
horas_noche = np.linspace(20, 28, 400)
reloj = (horas_noche % 24).reshape(-1, 1)

fig, ax = plt.subplots(figsize=(9.5, 4.2))
ax.plot(horas_noche, flujo_verdadero(reloj.ravel()), color=COLOR_VERDAD, lw=1.5, ls="--", label="Relación verdadera")
ax.plot(horas_noche, spline_normal.predict(reloj), color=COLOR_POLI, lw=2.2, label="Spline cúbico")
ax.plot(horas_noche, spline_periodico.predict(reloj), color=COLOR_SPLINE, lw=2.2, label="Spline cúbico periódico")
ax.axvline(24, color=COLOR_KNOT, ls=":", lw=1.4)
ax.set_xticks([20, 22, 24, 26, 28])
ax.set_xticklabels(["20:00", "22:00", "00:00", "02:00", "04:00"])
ax.set_xlabel("Hora (de la noche a la madrugada siguiente)")
ax.set_ylabel("Flujo (vehículos / min)")
ax.set_title("Las 24:00 y las 0:00 son el mismo instante: solo el periódico lo respeta")
ax.grid(True, alpha=0.25)
ax.legend(frameon=False)
plt.tight_layout()
plt.show()

pd.DataFrame(filas_periodico).round(3)

	Modelo	Columnas	RMSE CV	Flujo a las 24:00	Flujo a las 0:00	Diferencia
0	Spline cúbico	26	6.594	24.868	18.294	6.574
1	Spline cúbico periódico	23	6.537	20.070	20.070	0.000

El spline ordinario predice 24.9 vehículos/min a las 24:00 y 18.3 a las 0:00: una diferencia de 6.6. El periódico predice 20.1 en ambos extremos y une también las dos primeras derivadas. Esa igualdad describe un patrón diario periódico idealizado; dos fechas consecutivas reales no tienen por qué compartir el mismo tráfico.

El RMSE CV es 6.594 para el ordinario y 6.537 para el periódico. La comparación cambia simultáneamente frontera, ubicación de knots y número de columnas (26 frente a 23); por eso no aísla el efecto de la periodicidad. No optimizas otra rejilla para el periódico ni vuelves a abrir la prueba para elegirlo.

La restricción circular aporta coherencia cuando modelas un ciclo de 24 horas. Su conveniencia predictiva debe estudiarse con el diseño temporal que corresponda al uso final.

20.17 Preguntas de comprensión

1. Si la recta explica apenas el 6% de la variación del tráfico, ¿significa que la hora no importa? ¿El \(R^2\) está cerca de 0 por ruido o por mala especificación?
2. El polinomio de grado 12 tiene el RMSE in-sample más bajo de los polinomios (11.6). ¿Lo usarías para predecir las 23:30 de un día nuevo? ¿Qué hace entre las 0 y las 4, donde la curva verdadera es casi plana?
3. Con 15 knots internos, el RMSE in-sample es el más bajo (6.264). ¿Basta ese número para confiar en el modelo, o hace falta medirlo sobre datos que no vio?
4. El RMSE de validación cruzada del spline elegido (6.594) participó en la elección de n_knots=25. ¿Es una estimación honesta de su error en datos nuevos, y qué papel cumple la prueba reservada al principio?
5. ¿Por qué el capítulo usa dos muestras simuladas separadas, y es la regla de un error estándar una prueba de significancia?

TipOrientaciones

1. No. La hora determina la curva generadora; el \(R^2\) de 0.062 dice que una recta no puede representar los dos picos, y los residuos lo confirman, porque son positivos en los picos y negativos en la madrugada, el mediodía y la noche, con signo que cambia donde cambia el patrón. Es mala especificación, no ruido.
2. No. Reducir el RMSE in-sample no garantiza predecir bien en datos nuevos, y un polinomio global modifica todas las regiones a la vez: el grado 12 arranca en 38 vehículos por minuto a las 0:34, cuando la curva verdadera está cerca de 18, y el grado 5 se desploma hasta −45.6 a las 23:25.
3. No basta. El error de entrenamiento casi siempre baja al agregar knots, y 15 knots ganarían aunque estuvieran sobreajustando; la frontera la marca la validación cruzada, que en la segunda muestra encontró el mínimo en n_knots=25 y un repunte en 40 (RMSE de 7.713).
4. Puede ser optimista, porque ese mismo número sirvió para elegir. La prueba reservada de 45 observaciones evalúa el procedimiento ya fijado una sola vez, y dio un RMSE de 6.306 para el spline, frente a 10.485 del polinomio de grado 14 y 24.897 de la recta.
5. La primera muestra hace visibles las formas y sus fallos, y la segunda separa la selección y el reporte de esa exploración, con su prueba reservada antes de verla. La regla de un error estándar es una heurística que favorece la sencillez (aquí habría elegido n_knots=20); los folds comparten entrenamiento, y su error estándar ingenuo no tiene garantía general.

20.18 Referencias para profundizar

Splines y suavización

• de Boor, C. (2001). A Practical Guide to Splines (ed. revisada). Springer. La referencia clásica sobre B-splines.
• Hastie, T., Tibshirani, R. y Friedman, J. (2009). The Elements of Statistical Learning (2.ª ed.). Springer. Capítulo 5, expansiones de base y regularización.
• James, G., Witten, D., Hastie, T., Tibshirani, R. y Taylor, J. (2023). An Introduction to Statistical Learning with Applications in Python. Springer. Capítulo 7, más allá de la linealidad.
• Eilers, P. H. C. y Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science, 11(2), 89–102 (con discusión y respuesta hasta la página 121).
• Runge, C. (1901). Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten. Zeitschrift für Mathematik und Physik, 46, 224–243. El artículo original del fenómeno de Runge.

Usos en inteligencia artificial y analítica de datos

• Hastie, T. y Tibshirani, R. (1990). Generalized Additive Models. Chapman & Hall.
• Friedman, J. H. (1991). Multivariate adaptive regression splines. The Annals of Statistics, 19(1), 1–67.
• Balestriero, R. y Baraniuk, R. (2018). A spline theory of deep networks. Proceedings of the 35th International Conference on Machine Learning.

Documentación de las librerías

• scikit-learn, SplineTransformer: https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.SplineTransformer.html
• scikit-learn, ejemplo de variables de tiempo con splines periódicos: https://scikit-learn.org/stable/auto_examples/applications/plot_cyclical_feature_engineering.html

20.19 Lo que debes llevarte de este capítulo

Un spline da flexibilidad local con empalmes suaves, y cuánta flexibilidad darle se decide con validación cruzada, no a simple vista.

• Una recta obliga a un efecto constante: en el tráfico dejó \(R^2=0.062\) y residuos que repetían la forma de los picos.
• Un polinomio global modifica toda la curva a la vez: subir el grado bajó el RMSE in-sample hasta 11.6 y produjo flujos imposibles en los extremos, como −45.6 vehículos por minuto con grado 5; el ejemplo de Runge mostró oscilaciones de hasta −3.6 en una función que no baja de 0.04.
• Los tramos independientes saltan en los knots (de ~29 a ~120 a las 17:00), y un spline cúbico impone continuidad de \(f\), \(f'\) y \(f''\).
• Más knots dan más flexibilidad: el RMSE in-sample bajó de 18.923 con 3 knots internos a 6.264 con 15, sin decir nada sobre la generalización.
• La validación cruzada eligió n_knots=25 (RMSE de 6.594), y la prueba reservada confirmó el orden: 6.306 el spline, 10.485 el polinomio y 24.897 la recta.
• SplineTransformer con LinearRegression sigue siendo una regresión lineal sobre variables nuevas, y una hora circular pide un spline periódico, que unió las 0:00 y las 24:00 en 20.1 vehículos por minuto.

20.20 Puente hacia los siguientes capítulos

Los splines construyen una función flexible con polinomios locales y knots. El capítulo de métodos de kernel estima el comportamiento de una función usando directamente la información de los puntos cercanos a cada valor de \(X\), sin construir una base polinómica: una ventana alrededor de cada \(x\), un peso que decae con la distancia y una estimación local. En lugar de knots que colocar hay un ancho de ventana que elegir, y la misma validación cruzada lo calibra. Más adelante, el capítulo de modelos aditivos suma un spline por cada variable para modelar varias relaciones no lineales a la vez.