Generalizar: sobreajuste, validación y regularización

Tema 7 · Fase 2, Aprender de ejemplos · Toda la ruta

Cuaderno del tema. La teoría y los ejercicios con todo su código, listos para ejecutar en Colab sin instalar nada.

Los 15 bloques de código de la teoría, en el mismo orden. Si prefieres armar tu propio cuaderno, cópialos y córrelos en orden; en la ficha Teoría está la explicación de cada uno.

Datos de una regla conocida

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

# Semilla fija: cada vez que corras esto saldrán los mismos datos
rng = np.random.default_rng(1)
n = 40
x = rng.uniform(-3, 3, size=n)
# La regla «secreta»: una onda, más un ruido que nadie puede predecir
y = np.sin(x) + rng.normal(0, 0.3, size=n)
X = x.reshape(-1, 1)                  # scikit-learn quiere una tabla

# Primero se aparta la prueba (10 puntos): no se toca hasta el final
X_resto, X_prueba, y_resto, y_prueba = train_test_split(
    X, y, test_size=10, random_state=0)
# Del resto, 10 puntos para validación y 20 para entrenar
X_ent, X_val, y_ent, y_val = train_test_split(
    X_resto, y_resto, test_size=10, random_state=0)

print("Entrenamiento:", len(X_ent), "  validación:", len(X_val),
      "  prueba:", len(X_prueba))

Polinomios: una perilla de complejidad

from sklearn.preprocessing import PolynomialFeatures

# Las potencias de x, de 1 a 3, para x = 2 y x = 3
potencias = PolynomialFeatures(degree=3, include_bias=False)
print(potencias.fit_transform([[2], [3]]))
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error

def polinomio(grado):
    """Regresión polinomial: potencias de x, estandarizar, recta."""
    return make_pipeline(PolynomialFeatures(grado, include_bias=False),
                         StandardScaler(), LinearRegression())

# Una cuadrícula fina de x para dibujar las curvas
x_dib = np.linspace(-3, 3, 300).reshape(-1, 1)

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8),
                         sharex=True, sharey=True)
for eje, grado in zip(ejes.flat, [1, 3, 9, 15]):
    modelo = polinomio(grado).fit(X_ent, y_ent)
    e_ent = mean_squared_error(y_ent, modelo.predict(X_ent))
    e_val = mean_squared_error(y_val, modelo.predict(X_val))
    print(f"Grado {grado:2d}:  MSE entrenamiento = {e_ent:.3f}"
          f"   MSE validación = {e_val:.3f}")
    # Solo el grado que mejor generaliza va resaltado (color "C1")
    color = "C1" if grado == 3 else "C0"
    eje.scatter(X_ent, y_ent, s=14, color="C0")
    eje.scatter(X_val, y_val, s=14, facecolors="none", edgecolors="C2")
    eje.plot(x_dib, modelo.predict(x_dib), color=color)
    eje.set_title(f"Grado {grado}")
    eje.set_ylim(-2, 2)
for eje in ejes[1]:
    eje.set_xlabel("x")
for eje in ejes[:, 0]:
    eje.set_ylabel("y")
plt.show()

La curva en forma de U

grados = range(1, 16)
err_ent, err_val = [], []
for grado in grados:
    modelo = polinomio(grado).fit(X_ent, y_ent)
    err_ent.append(mean_squared_error(y_ent, modelo.predict(X_ent)))
    err_val.append(mean_squared_error(y_val, modelo.predict(X_val)))

mejor = grados[np.argmin(err_val)]
print("Grado con menor error de validación:", mejor)
for grado in [1, 3, 6, 12, 15]:
    print(f"Grado {grado:2d}:  entrenamiento {err_ent[grado - 1]:.3f}"
          f"   validación {err_val[grado - 1]:.3f}")

plt.plot(grados, err_ent, color="C2", marker="o")
plt.plot(grados, err_val, color="C0", marker="o")
plt.scatter([mejor], [min(err_val)], color="C1", s=160, zorder=3)
plt.yscale("log")
plt.ylim(0.006, 4)
# Nombres escritos junto a cada curva, en lugar de una leyenda
plt.text(15, 1.8, "Validación", ha="right")
plt.text(15, 0.008, "Entrenamiento", ha="right", color="C2")
plt.text(3.4, 0.016, "el mejor grado")
plt.xlabel("Grado del polinomio")
plt.ylabel("MSE (escala logarítmica)")
plt.title("Error según la complejidad")
plt.show()

Sesgo y varianza

def nuevos_datos(semilla, n=20):
    """Otra muestra de la misma regla: sen(x) más ruido."""
    r = np.random.default_rng(semilla)
    x_n = r.uniform(-3, 3, size=n)
    return x_n.reshape(-1, 1), np.sin(x_n) + r.normal(0, 0.3, size=n)

x0 = [[2.0]]                          # el punto donde miramos
for grado in [1, 3, 15]:
    preds = []
    for semilla in range(200):
        X_n, y_n = nuevos_datos(semilla)
        modelo = polinomio(grado).fit(X_n, y_n)
        preds.append(modelo.predict(x0)[0])
    preds = np.array(preds)
    print(f"Grado {grado:2d}:  mediana {np.median(preds):5.2f}"
          f"   la mitad central entre {np.percentile(preds, 25):5.2f}"
          f" y {np.percentile(preds, 75):5.2f}")

Validación cruzada k-fold

for semilla in range(5):
    Xe, Xv, ye, yv = train_test_split(X_resto, y_resto, test_size=10,
                                      random_state=semilla)
    modelo = polinomio(3).fit(Xe, ye)
    print(f"random_state={semilla}:  MSE de validación = "
          f"{mean_squared_error(yv, modelo.predict(Xv)):.3f}")
from sklearn.model_selection import KFold, cross_val_score

cinco = KFold(n_splits=5, shuffle=True, random_state=0)

resultados = []
for grado in grados:
    puntajes = cross_val_score(polinomio(grado), X_resto, y_resto,
                               cv=cinco, scoring="neg_mean_squared_error")
    resultados.append(-puntajes.mean())     # el MSE, con su signo

grado_cv = grados[np.argmin(resultados)]
for grado in [1, 2, 3, 4, 8, 15]:
    print(f"Grado {grado:2d}:  MSE de validación cruzada = "
          f"{resultados[grado - 1]:.3f}")
print("Grado elegido por validación cruzada:", grado_cv)
final = polinomio(grado_cv).fit(X_resto, y_resto)
error_prueba = mean_squared_error(y_prueba, final.predict(X_prueba))
print(f"Grado {grado_cv}, MSE con los datos de prueba: {error_prueba:.3f}")
# El MSE de la regla verdadera, sen(x), con estos mismos 10 puntos
print(f"MSE de sen(x) con los datos de prueba: "
      f"{mean_squared_error(y_prueba, np.sin(X_prueba[:, 0])):.3f}")

Regularización: castigar los pesos grandes

from sklearn.linear_model import Ridge

def polinomio_ridge(grado, alpha):
    """Como polinomio(), pero con Ridge en lugar de LinearRegression."""
    return make_pipeline(PolynomialFeatures(grado, include_bias=False),
                         StandardScaler(), Ridge(alpha=alpha))

for alpha in [0.001, 0.01, 0.1, 1, 10, 100]:
    modelo = polinomio_ridge(15, alpha).fit(X_ent, y_ent)
    e_ent = mean_squared_error(y_ent, modelo.predict(X_ent))
    e_val = mean_squared_error(y_val, modelo.predict(X_val))
    print(f"alpha = {alpha:7}:  entrenamiento {e_ent:.3f}"
          f"   validación {e_val:.3f}")

sin_reg = polinomio(15).fit(X_ent, y_ent)
con_reg = polinomio_ridge(15, 0.1).fit(X_ent, y_ent)
plt.scatter(X_ent, y_ent, s=18, color="C0")
plt.scatter(X_val, y_val, s=18, facecolors="none", edgecolors="C2")
plt.plot(x_dib, sin_reg.predict(x_dib), color="C2",
         label="Sin regularizar")
plt.plot(x_dib, con_reg.predict(x_dib), color="C1",
         label="Ridge, alpha = 0,1")
plt.ylim(-2, 2.6)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Grado 15, con y sin castigo")
# La leyenda va debajo del gráfico, para no tapar las curvas
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.25), ncol=2,
           frameon=False, handlelength=1.6, columnspacing=1.2)
plt.show()

Ridge y Lasso con datos reales

from sklearn.datasets import load_diabetes

diabetes = load_diabetes(as_frame=True)
nombres = list(diabetes.data.columns)
X_d = StandardScaler().fit_transform(diabetes.data)
y_d = diabetes.target.to_numpy()

alphas_r = np.logspace(-2, 4, 60)     # de 0,01 a 10 000
pesos_r = np.array([Ridge(alpha=a).fit(X_d, y_d).coef_ for a in alphas_r])

for i, nombre in enumerate(nombres):
    # Solo bmi va resaltado (color "C1"); las otras nueve, en "C2"
    if nombre == "bmi":
        plt.plot(alphas_r, pesos_r[:, i], color="C1", linewidth=3,
                 label="bmi")
    else:
        plt.plot(alphas_r, pesos_r[:, i], color="C2", linewidth=1.2)
plt.axhline(0, color="C0", linewidth=1)
plt.xscale("log")
plt.xlabel("alpha (escala logarítmica)")
plt.ylabel("Peso")
plt.title("Ridge: los pesos se encogen juntos")
plt.legend()
plt.show()

print("Pesos con alpha = 1000:",
      np.round(Ridge(alpha=1000).fit(X_d, y_d).coef_, 1))
from sklearn.linear_model import Lasso

alphas_l = np.logspace(-2, 2, 60)     # de 0,01 a 100
pesos_l = np.array([Lasso(alpha=a).fit(X_d, y_d).coef_ for a in alphas_l])

for i, nombre in enumerate(nombres):
    if nombre == "bmi":
        plt.plot(alphas_l, pesos_l[:, i], color="C1", linewidth=3,
                 label="bmi")
    else:
        plt.plot(alphas_l, pesos_l[:, i], color="C2", linewidth=1.2)
plt.axhline(0, color="C0", linewidth=1)
plt.xscale("log")
plt.xlabel("alpha (escala logarítmica)")
plt.ylabel("Peso")
plt.title("Lasso: los pesos se apagan uno a uno")
plt.legend()
plt.show()

for alpha in [0.1, 1, 5, 20]:
    lasso = Lasso(alpha=alpha).fit(X_d, y_d)
    vivos = [n for n, w in zip(nombres, lasso.coef_) if w != 0]
    print(f"alpha = {alpha:4}:  {len(vivos)} pesos distintos de cero:",
          ", ".join(vivos))

¿Por qué Lasso deja pesos en cero?

def perdida(w1, w2):
    """Un valle con el mínimo en (2, 1), más empinado en w1."""
    return 3 * (w1 - 2) ** 2 + (w2 - 1) ** 2

t = np.linspace(0, 2 * np.pi, 4000)
# El borde del círculo de radio 1 y el del rombo |w1| + |w2| = 1
circulo = np.cos(t), np.sin(t)
rombo = (np.cos(t) / (np.abs(np.cos(t)) + np.abs(np.sin(t))),
         np.sin(t) / (np.abs(np.cos(t)) + np.abs(np.sin(t))))

i_c = np.argmin(perdida(*circulo))
i_r = np.argmin(perdida(*rombo))
print(f"Ridge (círculo): w1 = {circulo[0][i_c]:.2f}"
      f"   w2 = {circulo[1][i_c]:.2f}")
print(f"Lasso (rombo):   w1 = {rombo[0][i_r]:.2f}"
      f"   w2 = {rombo[1][i_r]:.2f}")

W1, W2 = np.meshgrid(np.linspace(-1.5, 3, 200), np.linspace(-1.5, 2, 200))
# Una figura más alta, para que quepa la leyenda debajo
plt.figure(figsize=(4.8, 4.6))
plt.contour(W1, W2, perdida(W1, W2), levels=[0.5, 2, 4.5, 8, 12.5],
            colors="C2", linewidths=1)
plt.plot(*circulo, color="C0", linestyle="--", label="Ridge: círculo")
plt.plot(*rombo, color="C0", label="Lasso: rombo")
plt.scatter([circulo[0][i_c]], [circulo[1][i_c]], color="C0", s=80,
            zorder=3)
plt.scatter([rombo[0][i_r]], [rombo[1][i_r]], color="C1", s=140,
            zorder=4, label="Solución de Lasso: w2 = 0")
plt.scatter([2], [1], color="C0", marker="+", s=120)
plt.gca().set_aspect("equal")
plt.ylim(-1.5, 2)
plt.xlabel("w1")
plt.ylabel("w2")
plt.title("El rombo toca en una esquina")
# La leyenda va debajo del gráfico, para no tapar las elipses
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.3), fontsize=11,
           frameon=False, handlelength=1.6)
plt.show()

Curva de aprendizaje: ¿más datos ayudarían?

from sklearn.model_selection import learning_curve

tamanos, p_ent, p_val = learning_curve(
    make_pipeline(StandardScaler(), LinearRegression()),
    diabetes.data, y_d, train_sizes=[15, 20, 30, 50, 100, 200, 353],
    cv=cinco, scoring="neg_mean_squared_error",
    shuffle=True, random_state=0)
mse_ent = -p_ent.mean(axis=1)         # promedio de los cinco pliegues
mse_val = -p_val.mean(axis=1)

for t, e, v in zip(tamanos, mse_ent, mse_val):
    print(f"{t:3d} pacientes:  entrenamiento {e:6.0f}"
          f"   validación {v:6.0f}")

# Entrenamiento a rayas, para distinguirlo de validación
plt.plot(tamanos, mse_ent, color="C2", marker="o", linestyle="--",
         label="Entrenamiento")
plt.plot(tamanos, mse_val, color="C0", marker="o", label="Validación")
plt.xscale("log")
plt.yscale("log")
# Marcas legibles en el eje x: solo algunos tamaños
plt.xticks([15, 30, 100, 353], ["15", "30", "100", "353"])
plt.minorticks_off()
plt.xlabel("Pacientes para entrenar")
plt.ylabel("MSE (escala logarítmica)")
plt.title("Curva de aprendizaje: diabetes")
plt.legend()
plt.show()

La receta completa con diabetes

from sklearn.metrics import r2_score

Xd_ent, Xd_prueba, yd_ent, yd_prueba = train_test_split(
    diabetes.data, y_d, test_size=0.25, random_state=0)

# Como en el tema 6: entrenar y medir con los mismos pacientes
lineal = make_pipeline(StandardScaler(), LinearRegression())
lineal.fit(Xd_ent, yd_ent)
print(f"R² con los datos de entrenamiento: "
      f"{r2_score(yd_ent, lineal.predict(Xd_ent)):.3f}")

# Elegir alpha con validación cruzada, sin tocar la prueba
for alpha in [0.1, 1, 10, 100]:
    puntajes = cross_val_score(
        make_pipeline(StandardScaler(), Ridge(alpha=alpha)),
        Xd_ent, yd_ent, cv=cinco, scoring="r2")
    print(f"Ridge alpha = {alpha:5}:  R² de validación cruzada = "
          f"{puntajes.mean():.3f}")
# La prueba, una sola vez, con el modelo elegido
elegido = make_pipeline(StandardScaler(), Ridge(alpha=10))
elegido.fit(Xd_ent, yd_ent)
print(f"R² con los datos de prueba: "
      f"{r2_score(yd_prueba, elegido.predict(Xd_prueba)):.3f}")