Generalizar: sobreajuste, validación y regularización

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))

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()

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()

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}")

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}")

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")
# Más espacio arriba, para la leyenda
plt.ylim(-2, 3.2)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Grado 15, con y sin castigo")
plt.legend(loc="upper left")
plt.show()

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))

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))
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.8)
plt.xlabel("w1")
plt.ylabel("w2")
plt.title("El rombo toca en una esquina")
plt.legend(loc="upper left", fontsize=10, frameon=True)
plt.show()

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 folds
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}")

plt.plot(tamanos, mse_ent, color="C2", marker="o", 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()

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}")