Generalizar: sobreajuste, validación y regularización

1. Entrenamiento: es para lo único que existe.
2. Validación (o validación cruzada con los datos que no son de prueba): comparar modelos es una decisión.
3. Validación, por la misma razón: alpha es una perilla del modelo, como el grado.
4. Prueba, una sola vez y al final. Es el único grupo que no participó en ninguna decisión, así que es el único que da un número honesto.
5. Entrenamiento. Si la media y la desviación se calculan con todos los datos, información de la validación y de la prueba se cuela en el modelo. Por eso en la teoría el StandardScaler va dentro del pipeline: fit lo ajusta solo con los datos de entrenamiento.

• A, subajuste: los dos errores son altos y parecidos, muy por encima del ruido. El modelo no capta ni los datos que ve. Hace falta un modelo más flexible.
• B, sobreajuste fuerte: el error de entrenamiento es mucho menor que el ruido (está memorizando el ruido) y el de validación es enorme. Hay que simplificar el modelo o regularizar.
• C, buen ajuste: los dos errores están cerca del ruido y cerca entre sí.
• D, sobreajuste moderado: el error de entrenamiento parece bueno, pero el de validación es más del triple. Conviene regularizar un poco o probar un modelo más simple.

def elegir_grado(cv):
    """El grado de 1 a 15 con menor MSE de validación cruzada."""
    errores = []
    for grado in range(1, 16):
        puntajes = cross_val_score(polinomio(grado), X_resto, y_resto,
                                   cv=cv, scoring="neg_mean_squared_error")
        errores.append(-puntajes.mean())
    mejor = int(np.argmin(errores)) + 1
    return mejor, min(errores)

for nombre, cv in [("10 folds        ", KFold(10, shuffle=True,
                                             random_state=0)),
                   ("5 folds, otra mezcla", KFold(5, shuffle=True,
                                                  random_state=1))]:
    grado, error = elegir_grado(cv)
    print(f"{nombre}:  grado {grado}   MSE {error:.3f}")

10 folds        :  grado 3   MSE 0.042
5 folds, otra mezcla:  grado 3   MSE 0.048

Con 10 folds se entrenan 10 modelos por grado: 150 modelos en total.

Usó la prueba para elegir el grado. Al probar 15 modelos y quedarse con el que mejor le fue con esos 10 puntos, eligió en parte al que tuvo más suerte con ellos. El número que reporta ya no mide datos nuevos: es el mínimo de 15 intentos, y por eso es optimista.

X_nuevos, y_nuevos = nuevos_datos(999, n=2000)
reportado, real = [], []              # eligiendo con la prueba
reportado_3, real_3 = [], []          # grado 3 fijo, prueba una vez
for semilla in range(100):
    X_n, y_n = nuevos_datos(1000 + semilla, n=40)
    X_a, X_b = X_n[:30], X_n[30:]     # 30 para entrenar, 10 de prueba
    y_a, y_b = y_n[:30], y_n[30:]
    modelos = [polinomio(g).fit(X_a, y_a) for g in range(1, 16)]
    errores = [mean_squared_error(y_b, m.predict(X_b)) for m in modelos]
    elegido = modelos[int(np.argmin(errores))]
    reportado.append(min(errores))
    real.append(mean_squared_error(y_nuevos, elegido.predict(X_nuevos)))
    reportado_3.append(errores[2])    # el grado 3 es el tercer modelo
    real_3.append(
        mean_squared_error(y_nuevos, modelos[2].predict(X_nuevos)))

for nombre, rep, rea in [("Eligiendo con la prueba", reportado, real),
                         ("Grado 3 fijo           ", reportado_3, real_3)]:
    print(f"{nombre}:  reportado {np.mean(rep):.3f}"
          f"   real {np.mean(rea):.3f}"
          f"   (medianas {np.median(rep):.3f} y {np.median(rea):.3f})")

Eligiendo con la prueba:  reportado 0.091   real 1.201   (medianas 0.086 y 0.112)
Grado 3 fijo           :  reportado 0.105   real 0.111   (medianas 0.095 y 0.107)

Eligiendo con la prueba, el número reportado promete un error de alrededor de 0,09, pero el error real es mucho mayor: a veces el grado que tuvo suerte con 10 puntos es uno muy alto que se dispara con datos nuevos, y esos casos suben muchísimo el promedio. Incluso la mediana, que no se deja llevar por esos casos, sale más alta que lo reportado. Con el grado fijado de antemano, lo reportado y lo real coinciden bastante: la prueba usada una sola vez no engaña en promedio.

La forma correcta es la de la teoría: elegir con validación (o validación cruzada) sin tocar la prueba, y usar la prueba una sola vez, con el modelo ya elegido.

ridge = Ridge(alpha=1000).fit(X_d, y_d)
lasso = Lasso(alpha=5).fit(X_d, y_d)
for nombre, m in [("Ridge", ridge), ("Lasso", lasso)]:
    print(f"{nombre}:  pesos en cero = {np.sum(m.coef_ == 0)}"
          f"   el más chico = {np.min(np.abs(m.coef_)):.4f}")

Ridge:  pesos en cero = 0   el más chico = 0.0433
Lasso:  pesos en cero = 5   el más chico = 0.0000

a) Ridge, aun con un castigo enorme, no deja ningún peso en cero: el más chico es diminuto, pero no es cero. Lasso, con un castigo moderado, deja varios exactamente en cero.

b) La región de Ridge es un círculo, sin esquinas: la elipse del MSE lo toca en un punto cualquiera del borde, donde los pesos son chicos pero no nulos. La región de Lasso es un rombo, con esquinas sobre los ejes que sobresalen, y la elipse suele tocarlo primero en una esquina o en una arista, donde uno o varios pesos valen cero.

r2s = []
for semilla in range(10):
    Xa, Xb, ya, yb = train_test_split(diabetes.data, y_d, test_size=0.25,
                                      random_state=semilla)
    m = make_pipeline(StandardScaler(), Ridge(alpha=10)).fit(Xa, ya)
    r2s.append(r2_score(yb, m.predict(Xb)))
print("R² de prueba:", np.round(r2s, 2))
print(f"de {min(r2s):.2f} a {max(r2s):.2f}")

cv_todos = cross_val_score(make_pipeline(StandardScaler(), Ridge(alpha=10)),
                           diabetes.data, y_d, cv=cinco, scoring="r2")
print(f"Validación cruzada: {np.round(cv_todos, 2)}"
      f"   promedio {cv_todos.mean():.3f}")

R² de prueba: [0.36 0.44 0.45 0.43 0.47 0.53 0.45 0.51 0.47 0.6 ]
de 0.36 a 0.60
Validación cruzada: [0.33 0.46 0.53 0.51 0.61]   promedio 0.489

Con 111 pacientes de prueba, el R² cambia bastante según quiénes cayeron en la prueba: va de 0,36 a 0,60, y la semilla 0 fue justo la de peor suerte. Un solo número de prueba con pocos datos tiene un margen de error grande, y eso hay que decirlo al reportarlo. La validación cruzada promedia varias divisiones y da una idea de cuánto varía el resultado. Lo que no se vale es probar semillas hasta encontrar la que da el mejor número.

X_mal = X_d.copy()
X_mal[:, 2] = X_mal[:, 2] / 100         # bmi, en una unidad 100 veces mayor

for nombre, datos in [("estandarizado", X_d), ("bmi / 100    ", X_mal)]:
    m = Ridge(alpha=10).fit(datos, y_d)
    r2 = cross_val_score(Ridge(alpha=10), datos, y_d, cv=cinco,
                         scoring="r2").mean()
    # El efecto de bmi en las predicciones: peso por desviación de la columna
    efecto = m.coef_[2] * datos[:, 2].std()
    print(f"{nombre}:  efecto de bmi = {efecto:5.2f}   R² (CV) = {r2:.3f}")

estandarizado:  efecto de bmi = 24.60   R² (CV) = 0.489
bmi / 100    :  efecto de bmi =  0.08   R² (CV) = 0.421

Con bmi dividida por 100, para tener el mismo efecto el modelo necesitaría un peso 100 veces mayor, y Ridge castiga ese peso al cuadrado: le sale carísimo. El castigo prácticamente borra a bmi, la variable que más ayudaba, y el modelo empeora. Que una columna se castigue más que otra no debería depender de las unidades en que alguien la anotó: por eso se estandariza antes de regularizar.

from sklearn.datasets import fetch_california_housing

casas = fetch_california_housing(as_frame=True)
Xc_ent, Xc_prueba, yc_ent, yc_prueba = train_test_split(
    casas.data, casas.target, test_size=0.2, random_state=0)
casas_lineal = make_pipeline(StandardScaler(), LinearRegression())
casas_lineal.fit(Xc_ent, yc_ent)
print(f"R² entrenamiento: "
      f"{r2_score(yc_ent, casas_lineal.predict(Xc_ent)):.3f}")
print(f"R² prueba:        "
      f"{r2_score(yc_prueba, casas_lineal.predict(Xc_prueba)):.3f}")

R² entrenamiento: 0.609
R² prueba:        0.594

Los dos R² son casi iguales: con más de 16 000 zonas para entrenar y solo 8 pesos, una regresión lineal no tiene cómo memorizar. Su problema es el contrario, el sesgo: una función lineal de las 8 columnas es demasiado simple para los precios. Para mejorar haría falta un modelo más flexible (por ejemplo, los árboles del tema 11), no regularizar.

Xp_ent, Xp_prueba, yp_ent, yp_prueba = train_test_split(
    X_p, y_p, test_size=0.25, random_state=0)

comun = LinearRegression().fit(Xp_ent, yp_ent)
print(f"Regresión común:  R² entrenamiento "
      f"{comun.score(Xp_ent, yp_ent):.3f}"
      f"   R² prueba {comun.score(Xp_prueba, yp_prueba):.3f}")

alphas_p = np.logspace(-3, 0, 13)
r2_cv = [cross_val_score(Lasso(alpha=a), Xp_ent, yp_ent, cv=cinco,
                         scoring="r2").mean() for a in alphas_p]
alpha_p = alphas_p[np.argmax(r2_cv)]
print(f"Mejor alpha por validación cruzada: {alpha_p:.3f}"
      f"   (R² de validación cruzada {max(r2_cv):.3f})")

lasso_p = Lasso(alpha=alpha_p).fit(Xp_ent, yp_ent)
print(f"Lasso:  R² prueba {lasso_p.score(Xp_prueba, yp_prueba):.3f}")
elegidas = np.flatnonzero(lasso_p.coef_)
print("Columnas que usa Lasso:", elegidas)
print("Columnas que importan: ", np.sort(utiles))
colados = [j for j in elegidas if j not in utiles]
print("Pesos de las columnas que se colaron:",
      np.round(lasso_p.coef_[colados], 2))
print("Pesos de las 5 que importan (Lasso / verdad):")
for j in np.sort(utiles):
    print(f"  columna {j:2d}:  {lasso_p.coef_[j]:5.2f} / {w_verdad[j]:5.2f}")

Regresión común:  R² entrenamiento 0.951   R² prueba 0.811
Mejor alpha por validación cruzada: 0.178   (R² de validación cruzada 0.817)
Lasso:  R² prueba 0.880
Columnas que usa Lasso: [ 0  1  2  3  5  6  7 15 23 24 25 26]
Columnas que importan:  [ 3  5 24 25 26]
Pesos de las columnas que se colaron: [ 0.1   0.19 -0.24 -0.19  0.01  0.05 -0.04]
Pesos de las 5 que importan (Lasso / verdad):
  columna  3:  -1.34 / -1.76
  columna  5:   2.51 /  2.76
  columna 24:   1.55 /  2.06
  columna 25:   1.35 /  1.35
  columna 26:   2.40 /  2.62

• La regresión común sobreajusta: con 60 datos para entrenar y 30 pesos, explica el 95 % de la variación con los datos que vio, pero baja a 0,81 con la prueba.
• Lasso, con alpha elegido por validación cruzada, mejora la prueba a 0,88 y encuentra las 5 columnas que importan, con pesos cercanos a los verdaderos (casi todos algo más chicos: el castigo los encoge).
• Se le colaron 7 columnas inútiles, todas con pesos chicos (ninguno pasa de 0,25 en valor absoluto, contra más de 1,3 de las útiles). Con un alpha un poco mayor se apagarían, a costa de encoger más las útiles.
• Ridge encoge los pesos, pero los deja a los 30 distintos de cero: no te diría cuáles variables importan.