La validación cruzada estima qué tan bien generaliza un modelo apartando datos. El conjunto de datos se divide en un conjunto de entrenamiento y un conjunto de validación; el modelo aprende sobre el primero y se puntúa sobre el segundo. Repetir la división varias veces y promediar los scores da una estimación más confiable que una división única, y reduce el riesgo de elegir un modelo que solo memorizó los datos de entrenamiento. Vea los tutoriales de scikit-learn sobre validación cruzada.
De scikit-learn: una división única en conjuntos de entrenamiento y validación.
Un solo conjunto de validación es pequeño, así que su estimación del error es ruidosa. La validación cruzada sobre varios pliegues promedia la estimación sobre muchos subconjuntos apartados:
De scikit-learn.
En resumen, la validación cruzada cumple tres propósitos: estima el desempeño predictivo sobre datos no vistos, apoya la selección de modelos y el ajuste de hiperparámetros, y señala el sobreajuste.
Hay una trampa, y en las geociencias importa más que en casi cualquier otro campo. Nuestros datos están correlacionados en el tiempo y en el espacio: la posición GNSS de hoy es casi idéntica a la de ayer, y el ruido de una estación sísmica se parece al de su vecina. Cuando barajamos muestras correlacionadas al azar entre los pliegues, cada muestra de validación tiene gemelas cercanas sentadas en el conjunto de entrenamiento. La información se fuga del conjunto de validación hacia el entrenamiento, y los scores mienten. Esta lección demuestra la falla y su corrección sobre una serie sintética de desplazamiento GNSS con verdad de referencia (ground truth) conocida. Ediciones anteriores de este libro descargaban la estación P395 del Nevada Geodetic Laboratory — el mismo archivo global que sirve las estaciones de TLALOCNet en México —; la serie sintética conserva la misma física y añade una verdad de referencia contra la cual verificar. La segunda mitad de la lección pasa del tiempo al espacio y a los grupos: tablas multisitio donde la división debe respetar sitios, conglomerados y eventos.
1. Una serie de desplazamiento GNSS¶
El paquete mlgeo_synth genera una serie diaria de desplazamiento GNSS con componentes conocidas: una velocidad secular de 12 mm/año, carga estacional anual y semianual, un salto cosísmico de 25 mm en el día 1800 seguido de un decaimiento postsísmico logarítmico, y ruido coloreado (blanco + flicker + caminata aleatoria). Las componentes sin ruido regresan como columnas separadas, así que podemos comparar cualquier estimación con la verdad.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import mlgeo_synth
%matplotlib inlineg = mlgeo_synth.gnss_series(n_years=10.0, velocity_mm_yr=12.0, eq_day=1800, seed=42)
print(g.shape)
g.head()(3652, 5)
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
axes[0].plot(g["date"], g["disp_mm"], lw=0.5, color="0.3")
axes[0].set_ylabel("displacement (mm)")
axes[0].set_title("Observed displacement (with colored noise)")
axes[1].plot(g["date"], g["trend_mm"], label="trend (12 mm/yr)")
axes[1].plot(g["date"], g["seasonal_mm"], label="seasonal")
axes[1].plot(g["date"], g["eq_mm"], label="coseismic + postseismic")
axes[1].set_ylabel("displacement (mm)")
axes[1].set_title("Noise-free components")
axes[1].legend(loc="upper left")
axes[1].set_xlabel("date")
plt.tight_layout()
La serie observada es la suma de las tres componentes más el ruido. La tendencia registra el movimiento de las placas. El término estacional proviene de la carga hidrológica y atmosférica. El término del sismo es un salto más un decaimiento lento a medida que la falla se relaja. El ruido no es blanco: las componentes flicker y de caminata aleatoria hacen que los errores de días cercanos estén fuertemente correlacionados, que es exactamente lo que más adelante romperá la validación cruzada aleatoria.
Calentamiento: ajustar la velocidad secular. Antes de cualquier validación cruzada, ajustamos una línea recta a toda la serie, como haría un geodesta para estimar la velocidad de placa.
from sklearn.linear_model import LinearRegression
t_years = ((g["date"] - g["date"].iloc[0]).dt.days / 365.25).to_numpy().reshape(-1, 1)
d = g["disp_mm"].to_numpy()
reg = LinearRegression()
reg.fit(t_years, d)
d_fit = reg.predict(t_years)
print(f"Fitted velocity: {reg.coef_[0]:.2f} mm/yr")
print("True velocity: 12.00 mm/yr")
plt.figure(figsize=(10, 3.5))
plt.plot(g["date"], d, lw=0.5, color="0.5", label="data")
plt.plot(g["date"], d_fit, color="C3", lw=2, label="linear fit")
plt.ylabel("displacement (mm)")
plt.legend()
plt.title("Linear velocity fit on the full series")Fitted velocity: 19.52 mm/yr
True velocity: 12.00 mm/yr
La pendiente ajustada no da los 12 mm/año verdaderos. El salto cosísmico y el decaimiento postsísmico empujan la línea hacia arriba, y el ruido coloreado añade una deriva de período largo que una línea recta absorbe en parte. Un modelo puede ser simple, estar bien ajustado y aun así estar sesgado cuando la física de los datos es más rica que el modelo.
2. Una tarea supervisada de pronóstico¶
Las preguntas de validación cruzada se agudizan cuando el modelo hace algo operativo. Aquí va una: dada la historia reciente de la estación, ¿dónde estará el mes próximo? Para cada día ancla construimos características a partir de los últimos 60 días y predecimos el desplazamiento medio de los 30 días siguientes.
Características por día ancla:
- media de los últimos 60 días
- media de los últimos 10 días
- último valor observado
- pendiente lineal sobre los últimos 60 días (
np.polyfit) - seno y coseno del día del año (para codificar la estación del año)
Objetivo: desplazamiento medio sobre los 30 días siguientes.
Deslizamos el día ancla hacia adelante con un paso de 5 días, lo que da unas 700 muestras. Note lo que hace esta construcción: días ancla vecinos comparten casi todos los datos de su ventana. La muestra i y la muestra i+1 se solapan en 55 de los 60 días de historia y en 25 de los 30 días del objetivo. Las muestras están fuertemente correlacionadas por construcción, tal como lo estarían para cualquier serie temporal geofísica en ventanas.
disp = g["disp_mm"].to_numpy()
doy = g["date"].dt.dayofyear.to_numpy()
rows = []
for i in range(60, len(disp) - 30, 5): # need 60 days history, 30 days future
hist60 = disp[i - 60:i]
rows.append({
"mean_60d": hist60.mean(),
"mean_10d": disp[i - 10:i].mean(),
"last_val": disp[i - 1],
"slope_60d": np.polyfit(np.arange(60), hist60, 1)[0],
"doy_sin": np.sin(2 * np.pi * doy[i] / 365.25),
"doy_cos": np.cos(2 * np.pi * doy[i] / 365.25),
"target": disp[i:i + 30].mean(),
})
feat = pd.DataFrame(rows)
X = feat.drop(columns="target").to_numpy()
y = feat["target"].to_numpy()
print(f"{len(feat)} samples, {X.shape[1]} features")
feat.head()713 samples, 6 features
El modelo de todos los experimentos siguientes es el mismo bosque aleatorio, de modo que cualquier diferencia en los scores proviene de la división, no del modelo.
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor(n_estimators=200, random_state=42)3. La mentira optimista: la validación cruzada aleatoria¶
Primero, la receta estándar que enseña la mayoría de los tutoriales: barajar las muestras en 5 pliegues y puntuar con cross_val_score. Probamos dos esquemas barajados, KFold(shuffle=True) y ShuffleSplit, y reportamos R² y el error absoluto medio.
from sklearn.model_selection import cross_val_score, KFold, ShuffleSplit
def cv_scores(model, X, y, cv, name):
r2 = cross_val_score(model, X, y, cv=cv, scoring="r2")
mae = -cross_val_score(model, X, y, cv=cv, scoring="neg_mean_absolute_error")
print(f"{name:22s} R2 = {r2.mean():6.3f} +/- {r2.std():.3f} MAE = {mae.mean():5.2f} mm")
return r2.mean(), mae.mean()
kf_shuffled = KFold(n_splits=5, shuffle=True, random_state=42)
ss = ShuffleSplit(n_splits=5, test_size=0.2, random_state=42)
results = {}
results["KFold (shuffled)"] = cv_scores(model, X, y, kf_shuffled, "KFold (shuffled)")
results["ShuffleSplit"] = cv_scores(model, X, y, ss, "ShuffleSplit")KFold (shuffled) R2 = 1.000 +/- 0.000 MAE = 0.61 mm
ShuffleSplit R2 = 1.000 +/- 0.000 MAE = 0.58 mm
Los scores se ven excelentes: R² cercano a 1, errores de uno o dos milímetros. Si esto fuera un artículo, el modelo se declararía un éxito. No lo es. La sección siguiente muestra cuánto puntúa el mismo modelo cuando la división respeta el tiempo.
4. Divisiones honestas para series temporales¶
TimeSeriesSplit entrena sobre el pasado y valida sobre el futuro, que es como un modelo de pronóstico se usa en realidad. Como caso intermedio probamos también KFold sin barajar, que valida sobre bloques contiguos de tiempo: sin duplicados cercanos a través de la frontera, salvo en los bordes de los bloques.
from sklearn.model_selection import TimeSeriesSplit
kf_blocks = KFold(n_splits=5, shuffle=False)
tss = TimeSeriesSplit(n_splits=5)
results["KFold (blocks)"] = cv_scores(model, X, y, kf_blocks, "KFold (blocks)")
results["TimeSeriesSplit"] = cv_scores(model, X, y, tss, "TimeSeriesSplit")
summary = pd.DataFrame(results, index=["mean R2", "mean MAE (mm)"]).T.round(3)
summaryKFold (blocks) R2 = -1.739 +/- 2.132 MAE = 11.56 mm
TimeSeriesSplit R2 = -4.271 +/- 3.709 MAE = 16.45 mm
El mismo modelo, los mismos datos, y los esquemas honestos reportan un skill (habilidad de pronóstico) mucho peor. Para ver por qué, graficamos qué muestras caen en entrenamiento y cuáles en validación para cada esquema. El índice de muestra corre de izquierda a derecha, que aquí es también el tiempo.
schemes = {"ShuffleSplit": ss, "KFold (shuffled)": kf_shuffled, "TimeSeriesSplit": tss}
fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
for ax, (name, cv) in zip(axes, schemes.items()):
for fold, (tr_idx, va_idx) in enumerate(cv.split(X)):
ax.scatter(tr_idx, np.full(len(tr_idx), fold), marker="|", s=60,
color="C0", label="train" if fold == 0 else None)
ax.scatter(va_idx, np.full(len(va_idx), fold), marker="|", s=60,
color="C1", label="validation" if fold == 0 else None)
ax.set_yticks(range(5))
ax.set_ylabel("fold")
ax.set_title(name)
ax.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
axes[-1].set_xlabel("sample index (time order)")
plt.tight_layout()
En los dos paneles barajados, las muestras de validación naranjas están salpicadas entre las muestras de entrenamiento azules. Cada muestra de validación tiene vecinas inmediatas en el conjunto de entrenamiento, y esas vecinas comparten 55 de sus 60 días de historia. El bosque no necesita pronosticar nada; interpola a sus vecinas, y el score mide memorización. TimeSeriesSplit (panel inferior) siempre valida sobre datos posteriores a todo aquello con lo que se entrenó. El modelo debe extrapolar hacia un futuro que nunca ha visto, y el score refleja el skill que tendría en despliegue.
Los bloques contiguos (KFold sin barajar) quedan en el medio: la fuga de datos (data leakage) solo ocurre cerca de los bordes de los bloques, así que los scores caen la mayor parte del camino hacia el número honesto. Las divisiones por bloques son la idea correcta siempre que no hay un único eje temporal, por ejemplo bloques espaciales de estaciones.
5. Todo score necesita un modelo de referencia y una barra de error¶
Los scores honestos de la sección 4 se ven mal. ¿Qué tan mal? Dos hábitos ponen cualquier score en contexto: compararlo con un modelo de referencia (baseline) que no requiere aprendizaje, y ponerle una barra de error.
El modelo de referencia natural para un pronóstico es la persistencia: el futuro se parece al presente. Aquí el pronóstico de persistencia de la media de los próximos 30 días es simplemente el último valor observado, una característica que ya calculamos (last_val). No cuesta nada y codifica la propiedad más fuerte de la serie: cambia lentamente.
from sklearn.metrics import mean_absolute_error, r2_score
persist = feat["last_val"].to_numpy()
rows_bl = []
abs_err = [] # per-sample validation errors in time order, for the bootstrap below
for fold, (tr, va) in enumerate(tss.split(X)):
fitted = RandomForestRegressor(n_estimators=200, random_state=42).fit(X[tr], y[tr])
pred = fitted.predict(X[va])
abs_err.append(np.abs(y[va] - pred))
rows_bl.append({
"model MAE (mm)": mean_absolute_error(y[va], pred),
"persistence MAE (mm)": mean_absolute_error(y[va], persist[va]),
"model R2": r2_score(y[va], pred),
"persistence R2": r2_score(y[va], persist[va]),
})
abs_err = np.concatenate(abs_err)
pd.DataFrame(rows_bl).rename_axis("fold").round(2)La comparación da una lección de humildad. La persistencia pronostica la media de 30 días con un error de unos 2 mm y puntúa un R² cercano a 0.9 en cada pliegue; el bosque es un orden de magnitud peor. La razón es la extrapolación: la estación sigue moviéndose a 12 mm/año, así que cada objetivo de validación queda por encima del rango de objetivos con los que el bosque entrenó, y un bosque no puede predecir fuera del rango de sus etiquetas de entrenamiento. La persistencia cabalga la tendencia gratis.
Note los valores negativos de R² del modelo. R² compara un modelo contra la constante que predice la media de los objetivos de validación, así que R² < 0 significa que el modelo es peor que predecir la media — ningún skill en absoluto. Cada vez que un artículo reporte skill de pronóstico sin un modelo de referencia de persistencia (o de climatología), pregunte cuánto habría puntuado ese modelo de referencia.
La barra de error. Un número solo, como «MAE = 16 mm», esconde qué tan incierta es la propia estimación. El bootstrap le pone un intervalo: remuestree los errores de validación con reemplazo, recalcule el MAE de cada remuestra y lea un intervalo de confianza de la dispersión. Para errores independientes basta con remuestrear errores individuales (el bootstrap simple). Nuestros errores están autocorrelacionados — los días ancla vecinos comparten la mayor parte de su ventana —, así que el remuestreo de errores individuales finge que hay más información independiente de la que hay, y el intervalo sale demasiado angosto. El bootstrap de bloques móviles remuestrea bloques contiguos (aquí 12 muestras, o 60 días con el paso de 5 días), de modo que cada remuestra conserva la correlación local.
rng = np.random.default_rng(0)
n = len(abs_err)
n_boot = 2000
# simple bootstrap: resample individual errors
simple = np.array([rng.choice(abs_err, n).mean() for _ in range(n_boot)])
# moving-block bootstrap: resample contiguous blocks of 12 samples (60 days)
L = 12
n_blocks = int(np.ceil(n / L))
starts = rng.integers(0, n - L + 1, size=(n_boot, n_blocks))
block = np.array([
np.concatenate([abs_err[s:s + L] for s in row])[:n].mean() for row in starts
])
print(f"pooled honest MAE: {abs_err.mean():.1f} mm")
print(f"simple bootstrap 95% CI: "
f"[{np.percentile(simple, 2.5):.1f}, {np.percentile(simple, 97.5):.1f}] mm")
print(f"moving-block bootstrap 95% CI: "
f"[{np.percentile(block, 2.5):.1f}, {np.percentile(block, 97.5):.1f}] mm")pooled honest MAE: 16.4 mm
simple bootstrap 95% CI: [15.5, 17.5] mm
moving-block bootstrap 95% CI: [13.2, 20.1] mm
El intervalo por bloques es aproximadamente el doble de ancho que el simple. Los mismos datos, el mismo estadístico; la única diferencia es el supuesto de independencia. Los errores correlacionados llevan menos muestras efectivas de las que su conteo sugiere, y el intervalo angosto del bootstrap simple es una forma más del optimismo del que trata toda esta lección. Reporte el intervalo que corresponda a la estructura de correlación de los errores, y repórtelo junto a la estimación puntual.
6. Validación cruzada dejando uno afuera (LOOCV)¶
La LOOCV (Leave-One-Out CV) divide los datos en conjuntos de entrenamiento y validación n veces, donde n es el número de puntos de datos. En cada ronda, el conjunto de entrenamiento es todo menos una muestra, y el conjunto de validación es esa única muestra apartada.
Ventajas: bajo sesgo respecto de los datos de entrenamiento, y el resultado es determinista, pues no hay división aleatoria que repetir. Desventaja: cuesta n ajustes del modelo, lo que resulta caro para todo lo que no sean conjuntos de datos pequeños y modelos baratos.
Lo demostramos sobre el problema simple de regresión de la velocidad, submuestreado a un día de cada 20 para que los n ajustes sigan siendo rápidos.
from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import mean_squared_error
t_sub = t_years[::20] # ~180 points
d_sub = d[::20]
print(f"{len(d_sub)} points -> {len(d_sub)} fits")
loo = LeaveOneOut()
vels, mse_val = [], []
for train_idx, val_idx in loo.split(t_sub):
t_train, t_val = t_sub[train_idx], t_sub[val_idx]
d_train, d_val = d_sub[train_idx], d_sub[val_idx]
reg = LinearRegression().fit(t_train, d_train)
vels.append(reg.coef_[0])
mse_val.append(mean_squared_error(d_val, reg.predict(t_val)))
vels, mse_val = np.array(vels), np.array(mse_val)
print(f"Velocity estimates: mean {vels.mean():.3f} mm/yr, std {vels.std():.4f} mm/yr")
print(f"Mean validation MSE: {mse_val.mean():.2f} mm^2")183 points -> 183 fits
Velocity estimates: mean 19.427 mm/yr, std 0.0146 mm/yr
Mean validation MSE: 106.93 mm^2
La velocidad apenas se mueve cuando se elimina un punto, así que la dispersión entre los n ajustes es diminuta. Dos advertencias. Primera, la LOOCV rara vez vale su costo: n ajustes para una estimación del error que k pliegues aproximan con 5 o 10. Segunda, la LOOCV no corrige la correlación. Las vecinas inmediatas de cada punto apartado están siempre en el conjunto de entrenamiento, así que para una serie autocorrelacionada la LOOCV está cerca de la división más optimista posible: es el problema de los pliegues barajados llevado a su límite.
7. Búsqueda de hiperparámetros bajo validación cruzada honesta¶
Los modelos tienen parámetros aprendidos de los datos (las divisiones de los árboles, los pesos de una regresión) e hiperparámetros fijados antes del entrenamiento (la profundidad del árbol, el tamaño de las hojas). El ajuste de hiperparámetros busca la configuración que puntúa mejor bajo validación cruzada, y es práctica estándar. Los enfoques usuales:
- Ajuste manual: ajustar a mano, guiado por el conocimiento del problema. Bueno para construir intuición; no es sistemático.
- Búsqueda en malla (grid search): evaluar cada combinación de una malla predefinida. Exhaustiva, pero cara a medida que la malla crece. En scikit-learn:
GridSearchCV. - Búsqueda aleatoria: extraer combinaciones de distribuciones predefinidas con un presupuesto fijo de iteraciones. Cubre espacios amplios con menos costo que una malla. En scikit-learn:
RandomizedSearchCV. - Optimización bayesiana: modelar el score como función de los hiperparámetros y gastar las evaluaciones donde la mejora parece probable (p. ej.,
optuna,scikit-optimize).
La elección del esquema de validación cruzada dentro de la búsqueda importa tanto como la búsqueda misma. Ajustar contra pliegues barajados selecciona el modelo que mejor memoriza a sus vecinas. Nosotros ajustamos contra TimeSeriesSplit, así que el ganador es el mejor pronosticador.
from sklearn.model_selection import GridSearchCV
param_grid = {"max_depth": [2, 4, 8, None], "min_samples_leaf": [1, 5, 20]}
grid = GridSearchCV(
RandomForestRegressor(n_estimators=100, random_state=42),
param_grid,
cv=TimeSeriesSplit(n_splits=5),
scoring="neg_mean_absolute_error",
)
grid.fit(X, y)
print("Best params:", grid.best_params_)
print(f"Best CV MAE: {-grid.best_score_:.2f} mm")Best params: {'max_depth': 4, 'min_samples_leaf': 5}
Best CV MAE: 15.99 mm
Para la búsqueda aleatoria, pase los objetos de distribución congelados de scipy.stats. Una versión anterior de este cuaderno premuestreaba con randint.rvs(size=10), lo que colapsa la búsqueda aleatoria sobre una lista fija de diez valores; pasar la distribución congelada permite que cada iteración extraiga valores frescos.
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint
distributions = {"max_depth": randint(1, 10), "min_samples_leaf": randint(1, 10)}
rand = RandomizedSearchCV(
RandomForestRegressor(n_estimators=100, random_state=42),
distributions,
n_iter=15,
random_state=0,
cv=TimeSeriesSplit(n_splits=5),
scoring="neg_mean_absolute_error",
)
rand.fit(X, y)
print("Best params:", rand.best_params_)
print(f"Best CV MAE: {-rand.best_score_:.2f} mm")
print(f"Grid search best MAE was {-grid.best_score_:.2f} mm")Best params: {'max_depth': 9, 'min_samples_leaf': 5}
Best CV MAE: 15.99 mm
Grid search best MAE was 15.99 mm
Las dos búsquedas aterrizan a una fracción de milímetro una de la otra: con solo dos hiperparámetros, 15 extracciones aleatorias sondean el espacio casi tan bien como una malla de 12 puntos, y la búsqueda aleatoria escala mejor cuando el espacio crece.
8. Espacio y grupos: la misma mentira sobre un mapa¶
El tiempo no es el único eje a lo largo del cual los datos geocientíficos están correlacionados. Las campañas de campo producen datos agrupados: varios sondeos en cada sitio, varios sitios a lo largo de cada camino o cuenca, y un objetivo que cabalga sobre un campo regional suave. Dos observaciones del mismo sitio son casi copias; dos sitios del mismo conglomerado son primos cercanos.
mlgeo_synth.multisite_table construye esa geometría con verdad de referencia: 6 conglomerados de 5 sitios cada uno, 10 observaciones repetidas por sitio (300 filas). El objetivo es una señal lineal en las columnas feat_* — la parte transportable del skill — más un campo regional suave (longitud de correlación de 25 km) evaluado en cada sitio, más ruido. El campo se devuelve como un invocable en truth, así que podemos dibujar el mapa que el modelo está memorizando en secreto.
from sklearn.model_selection import GroupKFold, StratifiedGroupKFold
sites, truth = mlgeo_synth.multisite_table(seed=13)
print(sites.shape)
sites.head()(300, 8)
gx = np.linspace(-5, 105, 220)
GX, GY = np.meshgrid(gx, gx)
F = truth["field"](GX.ravel(), GY.ravel()).reshape(GX.shape)
fig, ax = plt.subplots(figsize=(7.5, 6))
im = ax.pcolormesh(GX, GY, F, cmap="viridis", shading="auto")
fig.colorbar(im, ax=ax, label="regional field (target units)")
site_locs = sites.drop_duplicates("site_id")
for cl, grp in site_locs.groupby("cluster_id"):
ax.scatter(grp["x_km"], grp["y_km"], s=50, edgecolor="white",
linewidth=1.2, label=f"cluster {cl}")
ax.set_xlabel("x (km)")
ax.set_ylabel("y (km)")
ax.set_title("Regional field with clustered sites (10 observations per site)")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncols=6, title=None)
plt.tight_layout()
El mapa es el mecanismo. El campo varía suavemente sobre ~25 km, y los sitios se asientan en conglomerados apretados de unos pocos km de ancho, así que todos los sitios de un conglomerado comparten casi el mismo valor del campo — y las 10 observaciones de un sitio lo comparten exactamente. Una división barajada le entrega al modelo esas casi copias.
Ahora, la escalera de divisiones. El mismo bosque, los mismos datos, tres divisiones que apartan progresivamente más de la estructura compartida: KFold barajado, GroupKFold sobre site_id (sitios enteros apartados) y GroupKFold sobre cluster_id con 6 divisiones — dejar un conglomerado afuera, la versión discreta de una división espacial por bloques. Corremos la escalera dos veces: una sobre las columnas feat_* solas, y otra añadiendo las coordenadas x_km, y_km como características.
feat_cols = [c for c in sites.columns if c.startswith("feat_")]
y_site = sites["target"].to_numpy()
ladder = {
"KFold (shuffled)": (KFold(n_splits=5, shuffle=True, random_state=42), None),
"GroupKFold (site)": (GroupKFold(n_splits=5), sites["site_id"]),
"GroupKFold (cluster)": (GroupKFold(n_splits=6), sites["cluster_id"]),
}
rows_ladder = []
for feats_name, cols in [("features only", feat_cols),
("features + x,y", feat_cols + ["x_km", "y_km"])]:
X_site = sites[cols].to_numpy()
for rung, (cv, groups) in ladder.items():
r2 = cross_val_score(model, X_site, y_site, cv=cv, groups=groups, scoring="r2")
rows_ladder.append({"features": feats_name, "split": rung,
"mean R2": r2.mean(), "std R2": r2.std()})
pd.DataFrame(rows_ladder).round(3)Lea la escalera de arriba hacia abajo; cada peldaño hace una pregunta más difícil, y más honesta.
- KFold (barajado), R² ≈ 0.7. Cada observación apartada tiene hasta 9 hermanas de su propio sitio en el entrenamiento. El bosque reconoce el sitio por sus características y recuerda el valor del campo de ese sitio. El score responde: «¿qué tan bien puedo predecir otro sondeo en un sitio que ya muestreé?»
- GroupKFold por sitio, R² ≈ 0.35. Se apartan sitios enteros, así que la memorización del sitio desaparece — pero los sitios apartados todavía tienen vecinos de conglomerado en el entrenamiento que comparten el campo. El score responde: «¿un sitio nuevo dentro de un conglomerado ya muestreado?»
- Dejar-un-conglomerado-afuera, R² < 0. Se apartan conglomerados enteros y el modelo enfrenta una región que nunca ha visto. El valor del campo allí es incognoscible desde los datos de entrenamiento, y los desplazamientos del campo que el bosque memorizó lo desorientan activamente — peor que predecir la media (sección 5). El score responde la pregunta que un mapa regional de amenaza realmente plantea: «¿una región nueva?»
Añadir x_km, y_km agudiza el contraste. Bajo validación cruzada barajada el bosque ahora interpola el mapa casi a la perfección (R² ≈ 0.95) — una habilidad genuinamente útil dentro de la región muestreada. Bajo dejar-un-conglomerado-afuera, las mismas coordenadas empeoran las cosas: apuntan a partes del mapa que los datos de entrenamiento nunca restringieron. El mismo modelo, las mismas características, veredictos opuestos — porque interpolar entre sitios muestreados y extrapolar hacia una región nueva son afirmaciones distintas, y la división decide cuál afirmación respalda el score. En escenarios continuos sin conglomerados naturales, la misma idea se convierte en una división espacial con zona de amortiguamiento: excluya del entrenamiento todo lo que quede a menos de una longitud de correlación de los sitios de validación.
El cuaderno 4.5 (sección 4.6) cruza esta línea deliberadamente: un ensamble profundo entrenado sobre un rango composicional acotado se evalúa fuera de rango — incluido el único error que su desacuerdo no logra señalar.
9. Datos pequeños, conglomerados y desbalanceados: cuando la propia validación cruzada por grupos se rompe¶
Los datos agrupados de la práctica geotécnica y geológica suelen ser, además, pequeños y desbalanceados: unos pocos cientos de casos históricos, una clase positiva rara (licuefacción observada, deslizamiento ocurrido, mineralización presente) y positivos que se agrupan en el espacio porque el campo que los gobierna lo hace. Esa combinación es la norma, no el caso excepcional, dondequiera que las etiquetas provengan de casos históricos y no de levantamientos. multisite_table(binary=True) la reproduce: la misma geometría de sitios de 300 filas, con el objetivo latente umbralizado para que el 12 % de las filas sean positivas.
Agrupar por sitio sigue siendo obligatorio — pero ahora puede fallar en sus propios términos. GroupKFold reparte sitios enteros sin mirar las etiquetas, y como los positivos están agrupados en el espacio, un pliegue puede quedar sin ningún positivo. El AUC ROC no está definido sobre un pliegue con una sola clase.
import warnings
from sklearn.ensemble import RandomForestClassifier
cases, truth_b = mlgeo_synth.multisite_table(binary=True, seed=1)
X_case = cases[feat_cols].to_numpy()
y_case = cases["label"].to_numpy()
print(f"{y_case.sum()} positives in {len(y_case)} rows ({y_case.mean():.0%})")
clf = RandomForestClassifier(n_estimators=200, random_state=42)
gkf = GroupKFold(n_splits=5)
pos_per_fold = [int(y_case[va].sum())
for _, va in gkf.split(X_case, y_case, groups=cases["site_id"])]
print("positives per validation fold:", pos_per_fold)
with warnings.catch_warnings():
warnings.simplefilter("ignore") # sklearn warns about the undefined fold
auc_gkf = cross_val_score(clf, X_case, y_case, cv=gkf,
groups=cases["site_id"], scoring="roc_auc")
print("GroupKFold AUC per fold:", np.round(auc_gkf, 3))36 positives in 300 rows (12%)
positives per validation fold: [11, 5, 0, 18, 2]
GroupKFold AUC per fold: [0.834 0.98 nan 0.997 0.754]
El nan no es un error de scikit-learn; son los datos reportando que esta división por grupos dejó un pliegue con cero positivos, así que no hay curva ROC que calcular — y promediar los pliegues restantes cambia en silencio qué es lo que estima esa media. La reparación es StratifiedGroupKFold: mantener cada sitio intacto y balancear la fracción de positivos entre los pliegues.
sgkf = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=0)
pos_sgkf = [int(y_case[va].sum())
for _, va in sgkf.split(X_case, y_case, groups=cases["site_id"])]
print("positives per validation fold:", pos_sgkf)
auc_sgkf = cross_val_score(clf, X_case, y_case, cv=sgkf,
groups=cases["site_id"], scoring="roc_auc")
print(f"StratifiedGroupKFold AUC = {auc_sgkf.mean():.3f} +/- {auc_sgkf.std():.3f}")positives per validation fold: [8, 7, 7, 4, 10]
StratifiedGroupKFold AUC = 0.925 +/- 0.025
Ahora todos los pliegues tienen positivos y todos devuelven un número. Conserve en el reporte la dispersión entre pliegues: con 36 positivos repartidos en cinco partes, el AUC de cada pliegue descansa sobre un puñado de eventos, y la dispersión (o un intervalo bootstrap, sección 5) es tan informativa como la media.
10. Agrupar por evento: réplicas y explosiones de cantera¶
En sismología el grupo suele ser el evento. mlgeo_synth.event_station_table construye un conjunto de datos de juguete de movimiento del suelo: 60 sismos en 4 conglomerados (compactos en el espacio y en el tiempo, como las secuencias de sismo principal y réplicas), cada uno registrado por las mismas 15 estaciones — 900 registros. El logaritmo de la aceleración máxima del suelo sigue una relación de atenuación de juguete en magnitud y distancia, más un término por evento compartido por los 15 registros de un evento y un término por estación. El término por evento es la fuga de datos: una división aleatoria pone 12 registros de un evento en el entrenamiento y 3 en la validación, y el modelo recibe crédito por recordar el desplazamiento propio de ese evento.
gm, truth_gm = mlgeo_synth.event_station_table(seed=0)
X_gm = gm[["magnitude", "dist_km"]].to_numpy()
y_gm = gm["log10_pga"].to_numpy()
r2_rand = cross_val_score(model, X_gm, y_gm,
cv=KFold(n_splits=5, shuffle=True, random_state=42), scoring="r2")
r2_event = cross_val_score(model, X_gm, y_gm, cv=GroupKFold(n_splits=5),
groups=gm["event_id"], scoring="r2")
print(f"KFold (shuffled): R2 = {r2_rand.mean():.3f} +/- {r2_rand.std():.3f}")
print(f"GroupKFold (event_id): R2 = {r2_event.mean():.3f} +/- {r2_event.std():.3f}")KFold (shuffled): R2 = 0.845 +/- 0.014
GroupKFold (event_id): R2 = 0.600 +/- 0.150
Agrupar por evento cuesta cerca de un cuarto del skill aparente, y el número agrupado es el que predice el desempeño sobre el próximo sismo. La versión de clasificación de esta trampa es común: cientos de réplicas de un mismo sismo principal, o explosiones repetidas de una misma cantera, son casi copias unas de otras, y una división aleatoria le permite a un clasificador puntuar alto reconociendo la fuente en lugar del tipo de fuente. Esa es exactamente la fuga de datos que el recuadro de advertencia de 3.5 declara y tolera — las tablas de características de la tabla de clasificación del curso (leaderboard) no traen identificador de evento por el cual agrupar. Aquí los metadatos existen, así que el costo de esa concesión es medible: es la brecha entre las dos filas de arriba.
11. La elección de la división¶
Cada esquema de esta lección responde la misma pregunta planteada a una escala distinta de estructura compartida. Antes de confiar en cualquier score, pregúntese: ¿qué estructura comparten mis datos que la división debe respetar — tiempo, sitio, evento o espacio? Observaciones que comparten una ventana de tiempo piden TimeSeriesSplit o bloques contiguos; que comparten un sitio o un instrumento, GroupKFold por sitio; que comparten un evento fuente, GroupKFold por evento; que comparten una región, dejar-un-conglomerado-afuera o una división espacial con zona de amortiguamiento. La auditoría de fugas de datos de 2.13 hace esta pregunta antes de entrenar; la elección de la división es la misma auditoría aplicada a la evaluación. Y cuando faltan los metadatos necesarios para agrupar, dígalo y declare qué puede y qué no puede significar el score, como lo hace el recuadro de la tabla de clasificación del curso en 3.5.