Overfitting y el balance sesgo-varianza

Qué cubre este capítulo

Cada modelo de este libro tiene una perilla de flexibilidad. Un árbol de decisión tiene su profundidad, un polinomio tiene su grado, una red neuronal tiene su ancho. Sube la perilla y el modelo puede ajustarse a más formas; el error de entrenamiento baja, y sigue bajando mientras más le subas. La trampa es que el número que de verdad te importa — el error sobre datos que no has visto — no sigue bajando. Cae, toca fondo, y luego vuelve a subir. Esa curva en forma de U es la imagen más importante de todo el aprendizaje supervisado, y este capítulo se trata de dibujarla, entenderla y ponerle nombre a las dos fuerzas que la doblan.

El vehículo es la regresión polinomial, porque hace la perilla concreta: ajusta una recta, luego una parábola, luego una curva de grado 11, y observa cómo trabaja cada una sobre los mismos datos. Construimos todo desde cero en NumPy — la expansión de features, el ajuste por mínimos cuadrados, la métrica de error — y luego corremos una descomposición empírica de sesgo y varianza, reajustando sobre cientos de remuestreos bootstrap para medir cuánto del error viene de que el modelo es demasiado rígido para alcanzar la verdad (sesgo) y cuánto de que se sacude persiguiendo el ruido (varianza). La pieza central es una animación: barre el grado desde una línea recta hasta una ondulación salvaje, y mira cómo el error de entrenamiento cae hacia cero mientras el error de test dibuja su U. Este es el capítulo al que el resto del curso vuelve a apuntar una y otra vez — la regularización, la validación cruzada y cada hiperparámetro que llegues a ajustar son, todos, intentos de encontrar el fondo de esa U.

Un poco de historia

La tensión es vieja. En el siglo XIV, Guillermo de Ockham argumentó que entre explicaciones que compiten debes preferir la que tiene menos suposiciones — la navaja de Occam, el principio de que una hipótesis más simple tiene más probabilidades de ser correcta. Esa es toda la intuición detrás de combatir el overfitting, seiscientos años antes: un modelo que necesita cien ondulaciones para explicar tus datos probablemente está explicando el ruido, no la señal.

El planteamiento moderno es más joven y más afilado. La descomposición del error de predicción en sesgo y varianza fue folclore de la estimación estadística durante décadas, pero el artículo que la ancló al machine learning es "Neural Networks and the Bias/Variance Dilemma" de Stuart Geman, Élie Bienenstock y René Doursat, publicado en 1992 en Neural Computation. Su argumento era que un modelo flexible — una red neuronal grande, en su caso — puede llevar su sesgo casi a cero, pero paga esa flexibilidad en varianza: el ajuste oscila salvajemente dependiendo de qué muestra de entrenamiento le tocó ver. No puedes tener ambas cosas gratis a la vez, y elegir dónde sentarte en ese balance es, argumentaban, el problema central de aprender a partir de datos finitos. El nombre "dilema sesgo-varianza" viene de ese artículo, y treinta años después sigue siendo el marco al que todos recurren cuando un modelo se ve perfecto en el conjunto de entrenamiento y se desmorona en producción.

La intuición

Aquí está el planteamiento, reducido a una dimensión para que puedas verlo todo. Hay una función verdadera suave — una onda senoidal — y dibujé 40 muestras ruidosas a partir de ella. El modelo nunca ve la curva suave; solo ve los puntos. Su trabajo es recuperar la forma únicamente a partir de los puntos, y qué tanto se le permite intentarlo lo fija el grado del polinomio.

Tres ajustes, los mismos datos. Una recta de grado 1, el grado que el error de test va a elegir, y un polinomio de grado 11:

La línea naranja no tiene remedio. Una recta no puede doblarse, así que parte la diferencia y se pierde toda la forma — es demasiado rígida, y fallaría de la misma manera sin importar cuáles 40 puntos le entregaras. Esa rigidez es el sesgo. La curva morada es la falla opuesta: tiene libertad suficiente para perseguir puntos individuales, así que se sacude hacia arriba y abajo para pasar cerca de ruido que ni siquiera es parte de la señal verdadera. Dale otros 40 puntos y se sacudiría hacia otro lado completamente distinto. Esa sensibilidad es la varianza. La curva cian, grado 5, es la que encontró el seno sin perseguir los puntos — y el punto de este capítulo es que "grado 5" no es una corazonada: es el fondo medible de una curva que estamos a punto de dibujar.

La matemática

Fija una sola entrada de test $x$. El target ahí es la función verdadera más ruido: $y = f(x) + \varepsilon$, donde $\varepsilon$ tiene media cero y varianza $\sigma^2$ y es la aleatoriedad irreducible que ningún modelo puede predecir. Nuestro modelo ajustado $\hat{f}$ depende del conjunto de entrenamiento particular $D$ que nos tocó muestrear, así que es un objeto aleatorio; imagina muestrear muchos conjuntos de entrenamiento y reajustar cada vez. La cantidad que importa es el error cuadrático esperado en $x$, promediado sobre esa aleatoriedad. Se separa en exactamente tres piezas:

Lee los tres términos uno por uno. El sesgo es la brecha entre el ajuste promedio y la verdad: $\mathbb{E}[\hat{f}(x)]$ es la curva que obtendrías al promediar los ajustes de todos esos conjuntos de entrenamiento, y si ese promedio aún no alcanza a $f(x)$, el modelo es demasiado rígido para representar la verdad sin importar cuántos datos le des. La varianza es cuánto tiembla un ajuste individual alrededor de ese promedio de un conjunto de entrenamiento al siguiente — la dispersión de la nube de curvas. El ruido, $\sigma^2$, es el piso: incluso un modelo perfecto que clavara $f$ exactamente seguiría fallando por el ruido de cada punto nuevo, así que ninguna astucia lleva el error por debajo de él.

Todo el balance vive en los primeros dos términos, y jalan en direcciones opuestas conforme subes la flexibilidad. Un modelo rígido (grado bajo) tiene sesgo alto y varianza baja: no puede alcanzar la verdad, pero apenas se mueve cuando los datos cambian. Un modelo flexible (grado alto) tiene sesgo bajo y varianza alta: su ajuste promedio puede trazar cualquier forma, pero cualquier ajuste individual está a merced del ruido particular que vio. El error total es su suma más el piso, así que minimizarlo significa encontrar el grado donde su suma toca fondo — no donde alguno de los dos por sí solo es más pequeño. Esa suma es la U.

El alcance del marco, y sus límites

La descomposición sesgo-varianza es un lente, no un algoritmo, y su valor está en que convierte una preocupación vaga — "¿mi modelo está haciendo overfitting?" — en dos preguntas que puedes responder por separado. Si a tu modelo le va mal tanto en entrenamiento como en test, eso es sesgo: es demasiado simple, y necesitas un modelo más flexible o mejores features, no más datos. Si le va genial en entrenamiento y mal en test, eso es varianza: es demasiado flexible para los datos que tienes, y los arreglos son los opuestos — simplifica el modelo, regularízalo, o consigue más datos. Saber contra cuál de los dos estás peleando es la diferencia entre un arreglo que ayuda y uno que empeora las cosas, y yo recurro a este marco con casi cada modelo que no se está comportando.

Dos advertencias honestas. Primera: la separación limpia en tres partes es exacta solo para regresión con error cuadrático; para clasificación y otras funciones de pérdida la misma historia se sostiene en espíritu pero el álgebra no es tan ordenada, y "sesgo" y "varianza" se vuelven más cualitativos. Segunda — y esta es la interesante — la U estricta es una afirmación sobre el régimen clásico, donde el modelo tiene menos parámetros que datos. Las redes profundas modernas corren rutinariamente con muchos más parámetros que puntos de datos y, pasado el punto en que interpolan el conjunto de entrenamiento a la perfección, el error de test puede caer de nuevo — la curva de "double descent" que Belkin y sus coautores documentaron en 2019. Eso no deroga el balance; dice que la imagen más allá del umbral de interpolación es más rica que una sola U. Para todo lo de este libro, y para la mayoría de los modelos que vas a poner en producción, la U única es el modelo mental correcto, y es el que vale la pena internalizar primero.

Los datos

El dataset es sintético y con semilla fija, versionado como data/synthetic.csv, así que el capítulo entero se reproduce al bit. El proceso generativo es deliberadamente simple: muestrea la entrada $x$ uniformemente en $[0, 1]$, fija el target verdadero en $\sin(2\pi x)$ — un periodo completo, curvatura de verdad — y agrega ruido gaussiano con desviación estándar 0.25. El conjunto de entrenamiento son 40 de esos puntos; una muestra separada de 4,000 puntos del mismo proceso es el conjunto de test que califica la generalización. La varianza del ruido es 0.0625 en unidades del target al cuadrado, y ese número es el piso de la matemática de arriba: ningún modelo, por listo que sea, baja su error de test por debajo de él.

Elegí el seno a propósito. Una recta no puede seguirlo, así que el grado 1 está forzado a hacer underfitting y el sesgo es inconfundible. Un grado intermedio traza un periodo limpiamente. Un grado alto tiene más libertad de la que 40 puntos ruidosos pueden restringir, así que ondula — y como estoy ajustando sobre potencias crudas, la ondulación se pone genuinamente salvaje cerca de los bordes del intervalo, que es overfitting que puedes ver con tus propios ojos. Aquí están los datos por sí solos, los puntos y la curva verdadera que el modelo tiene que recuperar a partir de ellos:

Constrúyelo, una función a la vez

La regresión polinomial suena a modelo nuevo, pero es regresión lineal disfrazada. El truco es la expansión de features: un polinomio de grado $p$ en $x$ es un modelo lineal en las columnas $1, x, x^2, \ldots, x^p$. Construye esa matriz — la matriz de Vandermonde — y toda la maquinaria de mínimos cuadrados del capítulo de regresión lineal aplica sin cambios. El grado es simplemente cuántas columnas le entregas al solver, que es exactamente la perilla de flexibilidad:

def vandermonde(x, degree, x_lo=0.0, x_hi=1.0):
    """Expand a 1-D input into polynomial features: columns 1, u, u^2, ..., u^p.

    This is the whole trick that turns linear regression into curve fitting. A
    degree-p polynomial c0 + c1 u + c2 u^2 + ... + cp u^p is a *linear* model in
    the columns [1, u, u^2, ..., u^p], so if we build that design matrix we can
    fit a curve with the exact machinery we'd use to fit a line. Row i is one
    sample; column j is that sample raised to the power j. The result is
    (len(x), degree + 1) — degree sets how many columns, i.e. how flexible the
    model is allowed to be.

    One numerical nicety: we first map the input from its domain [x_lo, x_hi]
    onto [-1, 1]. Raw powers of x on [0, 1] get badly correlated at high degree
    (x^10 and x^11 are nearly identical there), which makes the least-squares
    system ill-conditioned; centering to [-1, 1] keeps the columns well separated
    so the fit stays honest up to the degrees we push it to. It changes the basis,
    not the model — the fitted curve is the same polynomial either way.
    """
    x = np.asarray(x, dtype=float)
    u = 2.0 * (x - x_lo) / (x_hi - x_lo) - 1.0
    return np.vander(u, N=degree + 1, increasing=True)

El único detalle es el reescalado a $[-1, 1]$ antes de tomar potencias. Las potencias crudas de un número en $[0, 1]$ se vuelven casi indistinguibles a grado alto — $x^{10}$ y $x^{11}$ son casi la misma columna — y eso hace que el sistema de mínimos cuadrados quede mal condicionado, lo cual es un problema numérico, no de modelado. Centrar la entrada mantiene las columnas separadas para que el ajuste se mantenga honesto hasta el grado 11. Cambia la base, no el modelo.

Con la matriz de diseño en mano, el ajuste es una sola resolución de mínimos cuadrados. Usamos np.linalg.lstsq, que minimiza $\lVert Vc - y \rVert^2$ vía una SVD — la ruta numéricamente estable, que importa precisamente porque las matrices de Vandermonde de grado alto son el caso mal condicionado:

def fit_poly(x, y, degree):
    """Least-squares fit of a degree-p polynomial. Returns the coefficient vector.

    Build the Vandermonde design matrix, then solve the linear least-squares
    problem min_c ||V c - y||^2 for the coefficients c. We use np.linalg.lstsq,
    which solves it through an SVD — the numerically stable way, because at high
    degree the Vandermonde columns get badly correlated (x^11 and x^12 look
    nearly identical on [0, 1]) and forming the normal equations directly would
    lose precision. Same objective as ordinary least squares; only the number of
    columns changes.
    """
    V = vandermonde(x, degree)
    coef, *_ = np.linalg.lstsq(V, y, rcond=None)
    return coef

La predicción reconstruye la misma expansión de features en los puntos de consulta y toma la suma ponderada. El grado queda implícito en la longitud del vector de coeficientes, así que quien llama nunca lo pasa dos veces:

def predict_poly(coef, x):
    """Evaluate the fitted polynomial at new inputs: yhat = V(x) @ coef.

    Rebuild the same feature expansion at the query points and take the weighted
    sum. degree is implied by len(coef) - 1, so a caller never has to pass it
    twice.
    """
    degree = len(coef) - 1
    return vandermonde(x, degree) @ coef

La métrica es el error cuadrático medio, y quiero ser preciso con su unidad porque todo el capítulo se cuenta en ella: el MSE es el promedio de los errores al cuadrado, así que su unidad es la unidad del target al cuadrado — "unidades del target²" en todo lo que sigue. Elevar al cuadrado castiga los errores grandes más fuerte que los pequeños y es lo que hace que el álgebra de la descomposición funcione:

def mse(y_true, y_pred):
    """Mean squared error, in target units squared.

    The average of the squared misses. Squaring keeps sign out of it and
    punishes big misses harder than small ones; the unit is the target's unit
    squared, which is what "target units^2" means everywhere in this chapter.
    This is the number the whole overfitting story is told in.
    """
    y_true = np.asarray(y_true, dtype=float)
    y_pred = np.asarray(y_pred, dtype=float)
    return float(np.mean((y_true - y_pred) ** 2))

Ahora la pieza que hace el trabajo real de este capítulo: la descomposición empírica de sesgo y varianza. La matemática de arriba es una esperanza sobre conjuntos de entrenamiento, y aquí estimamos esa esperanza de la manera honesta — remuestreando de verdad. Toma muchos remuestreos bootstrap del conjunto de entrenamiento, ajusta un polinomio de grado $p$ en cada uno, y predice sobre una malla de puntos de evaluación. La dispersión de esas predicciones es la varianza; la brecha entre su promedio y la función verdadera es el sesgo:

def bias_variance_decompose(x_train, y_train, x_eval, f_true, degree,
                            n_boot, noise_var, rng):
    """Empirically split test error into bias^2, variance, and noise at one degree.

    The recipe follows the decomposition literally. Draw n_boot bootstrap
    resamples of the training set (sample the rows with replacement), fit a
    degree-p polynomial on each, and predict at every point of a fixed evaluation
    grid x_eval. That gives a cloud of predictions at each grid point — one per
    resample.

      - bias^2  = mean over x of (average prediction - true f(x))^2.
                  How far the *typical* fit sits from the truth. High for a stiff
                  model that can't reach the target shape.
      - variance = mean over x of the variance of the predictions across resamples.
                  How much the fit jitters when the training data changes. High
                  for a flexible model that chases noise.
      - noise    = the irreducible term, the noise variance baked into the data.
                  No model can beat it; we pass it in because we know it here.

    Returns (bias2, variance, noise, total) where total = bias2 + variance +
    noise is the model's expected squared error on a fresh point.
    """
    n = len(x_train)
    preds = np.empty((n_boot, len(x_eval)))
    for b in range(n_boot):
        idx = rng.integers(0, n, size=n)           # bootstrap resample, with replacement
        coef = fit_poly(x_train[idx], y_train[idx], degree)
        preds[b] = predict_poly(coef, x_eval)

    mean_pred = preds.mean(axis=0)                  # average fit at each grid point
    bias2 = float(np.mean((mean_pred - f_true(x_eval)) ** 2))
    variance = float(np.mean(preds.var(axis=0)))
    noise = float(noise_var)
    total = bias2 + variance + noise
    return bias2, variance, noise, total

Esa función es el experimento sobre el que descansa todo el balance. Todo lo demás es graficar.

Míralo trabajar

Aquí está la imagen por la que existe este capítulo. Presiona play y el grado del polinomio barre de 1 a 11, un reajuste real por cuadro. El panel superior es la curva ajustada sobre los datos, con el seno verdadero punteado detrás. El panel inferior construye las dos curvas de error conforme el grado sube — el error de entrenamiento en cian, el error de test en morado — con un punto que marca el grado actual. El pie de imagen nombra el grado y ambos errores en unidades del target².

Mira los dos paneles juntos, porque la historia está en su desacuerdo. En grado 1 el ajuste es una recta que falla feo con el seno, y ambos errores son altos — eso es puro sesgo, el modelo demasiado rígido para doblarse. Conforme el grado sube al rango medio la curva se acomoda sobre el seno, y ambos errores caen juntos a su punto más bajo alrededor del grado 5: error de entrenamiento 0.041, error de test 0.061 en unidades del target². Entonces los paneles se separan. El error de entrenamiento sigue deslizándose — 0.038 en grado 11, un ajuste cada vez más pegado a los 40 puntos que puede ver. Pero el error de test da la vuelta y sube, suavemente al principio, luego con fuerza: 0.079 en grado 7, 0.154 en grado 8, y 0.817 en grado 11. Arriba puedes ver por qué — la curva deja de trazar el seno y empieza a lanzarse hacia puntos individuales, disparándose fuera del panel cerca de los bordes para pasar por puntos que son mayormente ruido. Esa brecha creciente entre las líneas cian y morada es el overfitting, dibujado en números reales. El modelo ya no está aprendiendo; está memorizando.

El panel inferior es una sola realización, eso sí — un conjunto de entrenamiento, un conjunto de test — y una sola realización es un poco dentada. La descomposición la suaviza en sus dos componentes. Reajusta sobre 400 remuestreos bootstrap en cada grado, mide el sesgo² y la varianza de las predicciones, suma el piso de ruido, y la U se resuelve en las dos fuerzas que tiene debajo:

Este es el dilema sesgo-varianza como gráfica, y el eje y es logarítmico porque los dos componentes abarcan órdenes de magnitud. El sesgo² (naranja) empieza alto — 0.123 en grado 1, la incapacidad de la recta rígida de alcanzar el seno — y colapsa en cuanto el modelo tiene grado suficiente para trazar la forma, aplanándose cerca de 0.003 desde el grado 3 en adelante. La varianza (morado) hace la imagen espejo: empieza casi en nada, 0.012 en grado 1, y sube conforme crece la flexibilidad, suavemente por los grados medios y luego explosivamente — 0.29 en grado 9, 1.6 en grado 11 — conforme el ajuste empieza a oscilar salvajemente entre remuestreos. Se cruzan alrededor del grado 4. El error total (cian) es su suma más el piso de ruido de 0.0625, y es la U: baja mientras el sesgo cae más rápido de lo que la varianza sube, sube en cuanto la varianza toma el control. El fondo de esa U, grado 3 aquí, es el punto dulce que todo este aparato existe para encontrar. A su izquierda estás limitado por sesgo; a su derecha, limitado por varianza.

La implementación completa

El archivo entero, sin librerías, de arriba a abajo — la expansión de features, el ajuste, la predicción, la métrica y la descomposición bootstrap. Este es el código que la animación y ambas gráficas corrieron de verdad:

"""Polynomial regression and the bias-variance decomposition, from scratch.

This chapter is about a single knob — model flexibility — and what it does to
generalization. The vehicle is polynomial regression: fit y with a polynomial of
degree p, turn p up, and watch the fit go from too stiff to draw the data (a
line through a curve) to so loose it threads every noisy point. Somewhere in the
middle is the fit that generalizes, and the whole point is that "somewhere in
the middle" is not a slogan — it's a measurable minimum of the test error.

The mechanism underneath is linear regression. A degree-p polynomial is a linear
model on an expanded feature set: instead of the single column x you fit on the
columns 1, x, x^2, ..., x^p. That expansion is a Vandermonde matrix, and once
you have it the fit is the same least-squares solve as ordinary regression. So
nothing here is new machinery — the flexibility all lives in how many columns
you hand the solver.

The second half is the bias-variance decomposition: refit the model on many
resampled training sets, look at the spread of the predictions, and split the
expected error into bias (how far the average fit sits from the truth), variance
(how much the fit jitters from sample to sample), and the irreducible noise. As
degree climbs, bias falls and variance rises; where they cross is the sweet spot.

Pure NumPy — no ML library in this file. Every function appears in the chapter
one region at a time via the include directives.
"""

import numpy as np


# region: vandermonde
def vandermonde(x, degree, x_lo=0.0, x_hi=1.0):
    """Expand a 1-D input into polynomial features: columns 1, u, u^2, ..., u^p.

    This is the whole trick that turns linear regression into curve fitting. A
    degree-p polynomial c0 + c1 u + c2 u^2 + ... + cp u^p is a *linear* model in
    the columns [1, u, u^2, ..., u^p], so if we build that design matrix we can
    fit a curve with the exact machinery we'd use to fit a line. Row i is one
    sample; column j is that sample raised to the power j. The result is
    (len(x), degree + 1) — degree sets how many columns, i.e. how flexible the
    model is allowed to be.

    One numerical nicety: we first map the input from its domain [x_lo, x_hi]
    onto [-1, 1]. Raw powers of x on [0, 1] get badly correlated at high degree
    (x^10 and x^11 are nearly identical there), which makes the least-squares
    system ill-conditioned; centering to [-1, 1] keeps the columns well separated
    so the fit stays honest up to the degrees we push it to. It changes the basis,
    not the model — the fitted curve is the same polynomial either way.
    """
    x = np.asarray(x, dtype=float)
    u = 2.0 * (x - x_lo) / (x_hi - x_lo) - 1.0
    return np.vander(u, N=degree + 1, increasing=True)
# endregion


# region: fit_poly
def fit_poly(x, y, degree):
    """Least-squares fit of a degree-p polynomial. Returns the coefficient vector.

    Build the Vandermonde design matrix, then solve the linear least-squares
    problem min_c ||V c - y||^2 for the coefficients c. We use np.linalg.lstsq,
    which solves it through an SVD — the numerically stable way, because at high
    degree the Vandermonde columns get badly correlated (x^11 and x^12 look
    nearly identical on [0, 1]) and forming the normal equations directly would
    lose precision. Same objective as ordinary least squares; only the number of
    columns changes.
    """
    V = vandermonde(x, degree)
    coef, *_ = np.linalg.lstsq(V, y, rcond=None)
    return coef
# endregion


# region: predict_poly
def predict_poly(coef, x):
    """Evaluate the fitted polynomial at new inputs: yhat = V(x) @ coef.

    Rebuild the same feature expansion at the query points and take the weighted
    sum. degree is implied by len(coef) - 1, so a caller never has to pass it
    twice.
    """
    degree = len(coef) - 1
    return vandermonde(x, degree) @ coef
# endregion


# region: mse
def mse(y_true, y_pred):
    """Mean squared error, in target units squared.

    The average of the squared misses. Squaring keeps sign out of it and
    punishes big misses harder than small ones; the unit is the target's unit
    squared, which is what "target units^2" means everywhere in this chapter.
    This is the number the whole overfitting story is told in.
    """
    y_true = np.asarray(y_true, dtype=float)
    y_pred = np.asarray(y_pred, dtype=float)
    return float(np.mean((y_true - y_pred) ** 2))
# endregion


# region: bias_variance
def bias_variance_decompose(x_train, y_train, x_eval, f_true, degree,
                            n_boot, noise_var, rng):
    """Empirically split test error into bias^2, variance, and noise at one degree.

    The recipe follows the decomposition literally. Draw n_boot bootstrap
    resamples of the training set (sample the rows with replacement), fit a
    degree-p polynomial on each, and predict at every point of a fixed evaluation
    grid x_eval. That gives a cloud of predictions at each grid point — one per
    resample.

      - bias^2  = mean over x of (average prediction - true f(x))^2.
                  How far the *typical* fit sits from the truth. High for a stiff
                  model that can't reach the target shape.
      - variance = mean over x of the variance of the predictions across resamples.
                  How much the fit jitters when the training data changes. High
                  for a flexible model that chases noise.
      - noise    = the irreducible term, the noise variance baked into the data.
                  No model can beat it; we pass it in because we know it here.

    Returns (bias2, variance, noise, total) where total = bias2 + variance +
    noise is the model's expected squared error on a fresh point.
    """
    n = len(x_train)
    preds = np.empty((n_boot, len(x_eval)))
    for b in range(n_boot):
        idx = rng.integers(0, n, size=n)           # bootstrap resample, with replacement
        coef = fit_poly(x_train[idx], y_train[idx], degree)
        preds[b] = predict_poly(coef, x_eval)

    mean_pred = preds.mean(axis=0)                  # average fit at each grid point
    bias2 = float(np.mean((mean_pred - f_true(x_eval)) ** 2))
    variance = float(np.mean(preds.var(axis=0)))
    noise = float(noise_var)
    total = bias2 + variance + noise
    return bias2, variance, noise, total
# endregion

La versión de librería

Nadie arma una matriz de Vandermonde a mano en producción. scikit-learn divide el mismo trabajo en dos pasos componibles: PolynomialFeatures construye la matriz de diseño expandida, y LinearRegression resuelve el ajuste de mínimos cuadrados sobre ella. Encadénalos en un pipeline y tienes exactamente nuestro fit_poly más predict_poly:

def sklearn_fit(x, y, degree):
    """Fit a degree-p polynomial with PolynomialFeatures + LinearRegression.

    Returns the fitted pipeline. include_bias=True gives the constant column, so
    LinearRegression runs with fit_intercept=False — the twin of our fit_poly,
    same design matrix, same least-squares solve. We rescale the input to [-1, 1]
    first, matching impl.vandermonde, so the coefficients line up one for one.
    """
    model = make_pipeline(
        PolynomialFeatures(degree=degree, include_bias=True),
        LinearRegression(fit_intercept=False),
    )
    model.fit(_rescale(x), y)
    return model

La predicción corre el mismo pipeline hacia adelante:

def sklearn_predict(model, x):
    """Predict with the fitted pipeline at new inputs."""
    return model.predict(_rescale(x))

La única sutileza es el intercepto. PolynomialFeatures ya emite la columna constante, así que ponemos fit_intercept=False para evitar un término de sesgo duplicado, y reescalamos la entrada a $[-1, 1]$ primero para empatar exactamente con nuestra Vandermonde. Con esas dos decisiones el vector de coeficientes del pipeline se alinea uno a uno con el nuestro, y como ambos usan el mismo solver de mínimos cuadrados vía SVD por debajo, aterrizan en la misma curva en cada grado — los coeficientes coinciden a alrededor de 4e-15 y las predicciones a 1e-12. No hay nada en lo que puedan estar en desacuerdo; es la misma álgebra lineal.

Desde cero contra la librería

Como las dos implementaciones son el mismo modelo, el duelo no es realmente desde cero contra librería — es una verificación de que construimos la referencia correctamente, y así fue. En el punto dulce, grado 5, ambas alcanzan un MSE de test de 0.061 en unidades del target², idéntico a la precisión citada arriba. La comparación interesante es de la que trata todo el capítulo: el costo de elegir el grado equivocado. Aquí está el error de test sobre datos apartados en el punto dulce junto al mismo error en el grado con overfitting, ambos en unidades del target²:

El modelo de grado 11 con overfitting anota 0.817 contra el 0.061 del punto dulce — 13 veces el error, sobre los mismos datos, de la misma familia de modelos, comprado puramente por girar la perilla de flexibilidad de más. Y aquí está la parte que hace al overfitting tan traicionero: sobre los datos de entrenamiento el modelo de grado 11 se ve mejor, 0.038 contra 0.041, así que si lo juzgaras por los números que te entrega sobre los datos con los que entrenó, mandarías a producción el peor modelo y te sentirías bien haciéndolo. La única defensa es calificar sobre datos que el modelo nunca vio, y por eso cada número de este capítulo que importa es un número de test, y por eso el conjunto de test es 100 veces más grande que el de entrenamiento — para que la U que traza sea lo bastante suave para confiar en ella. Las cifras exactas viven en results.json, regenerado cada vez que el código cambia, para que la prosa y las gráficas no puedan desviarse de lo que el código produjo.

Conclusiones

La curva en U es universal. No le importa que hayamos usado polinomios — cambia por la profundidad de un árbol, el ancho de una red, el número de vecinos en k-NN, el número de pasos de entrenamiento antes de detenerte, y obtienes la misma forma: el error de entrenamiento deslizándose hacia cero mientras el error de test cae, toca fondo y sube. Internaliza esa imagen y una cantidad enorme del machine learning práctico se vuelve una sola pregunta: ¿dónde está el fondo de mi U, y de qué lado de él estoy? El underfitting y el overfitting no son dos problemas sin relación; son los dos extremos de una misma perilla, y todo el juego es aterrizar en medio.

Una vez que sabes de qué lado estás, las jugadas son opuestas, y este es el premio de mantener sesgo y varianza separados en tu cabeza. Si estás limitado por sesgo — mal en entrenamiento y en test por igual — más datos no van a ayudar; necesitas un modelo más flexible o mejores features. Si estás limitado por varianza — genial en entrenamiento, pobre en test — los arreglos intercambian, todos, un poco de varianza por un poco de sesgo: simplifica el modelo, consigue más datos para sujetar al modelo flexible, o regularízalo, que es el siguiente hilo de este libro. Ridge y lasso son exactamente este intercambio hecho continuo — en vez de elegir un grado, agregas una penalización que jala los coeficientes hacia cero y ajustas su fuerza para deslizarte a lo largo de la U. La validación cruzada es cómo encuentras el fondo sin quemar un conjunto apartado. Todas son la misma idea con ropa distinta.

Una advertencia para llevar contigo, para que no te sorprenda después: la U única y ordenada es el régimen clásico, donde el modelo es más pequeño que los datos. Los modelos muy grandes — las redes profundas modernas en especial — pueden pasar el punto de ajustar el conjunto de entrenamiento a la perfección y ver el error de test mejorar de nuevo, el fenómeno de double descent. Es real y vale la pena conocerlo, pero no cambia el consejo para nada de lo que vayas a construir en este curso. Para polinomios, árboles, bosques, modelos lineales y prácticamente cualquier problema tabular que te vayas a encontrar, la U única es el mapa, y encontrar su fondo es el trabajo.