La validation croisée estime la capacité d’un modèle à généraliser en tenant des données à l’écart. Le jeu de données est découpé en un ensemble d’entraînement et un ensemble de validation ; le modèle apprend sur le premier et est noté sur le second. Répéter le découpage plusieurs fois et moyenner les scores donne une estimation plus fiable qu’un découpage unique, et réduit le risque de retenir un modèle qui n’a fait que mémoriser les données d’entraînement. Voir les tutoriels scikit-learn sur la validation croisée.
D’après scikit-learn : un découpage unique en ensembles d’entraînement et de validation.
Un ensemble de validation unique est petit, donc son estimation de l’erreur est bruitée. La validation croisée sur plusieurs plis moyenne l’estimation sur de nombreux sous-ensembles tenus à l’écart :
D’après scikit-learn.
En résumé, la validation croisée sert trois objectifs : elle estime la performance prédictive sur des données non vues, elle appuie la sélection de modèle et le réglage des hyperparamètres, et elle signale le surapprentissage.
Il y a un piège, et il compte plus en géosciences que presque partout ailleurs. Nos données sont corrélées dans le temps et dans l’espace : une position GNSS d’aujourd’hui est quasi identique à celle d’hier, et le bruit d’une station sismique ressemble à celui de sa voisine. Le signal hydrologique saisonnier des stations RENAG des Alpes, rencontré à la leçon 1.7, en est l’illustration : le déplacement d’un jour au suivant, et d’une station alpine à sa voisine, porte le même cycle de charge — deux échantillons voisins n’apportent pas deux informations indépendantes. Quand nous mélangeons au hasard des échantillons corrélés dans des plis, chaque échantillon de validation a des quasi-jumeaux assis dans l’ensemble d’entraînement. L’information fuit de l’ensemble de validation vers l’entraînement, et les scores mentent. Cette leçon démontre la défaillance et son remède sur une série de déplacement GNSS synthétique dont la vérité terrain est connue. Les éditions antérieures de ce livre téléchargeaient la station P395 du Nevada Geodetic Laboratory ; la série synthétique conserve la même physique et ajoute une vérité terrain à laquelle nous pouvons nous comparer. La seconde moitié de la leçon passe du temps à l’espace et aux groupes : des tableaux multi-sites où le découpage doit respecter les sites, les groupes spatiaux et les événements.
1. Une série de déplacement GNSS¶
Le paquet mlgeo_synth engendre une série journalière de déplacement GNSS aux composantes connues : une vitesse séculaire de 12 mm/an, une charge saisonnière annuelle et semi-annuelle, un saut cosismique de 25 mm au jour 1800 suivi d’une décroissance postsismique logarithmique, et un bruit coloré (blanc + scintillement + marche aléatoire). Les composantes sans bruit sont renvoyées dans des colonnes séparées, ce qui nous permet de comparer n’importe quelle estimation à la vérité.
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 série observée est la somme des trois composantes plus le bruit. La tendance enregistre le mouvement des plaques. Le terme saisonnier vient de la charge hydrologique et atmosphérique. Le terme sismique est un saut suivi d’une lente décroissance à mesure que la faille se relaxe. Le bruit n’est pas blanc : les composantes de scintillement et de marche aléatoire rendent les erreurs de jours voisins fortement corrélées — précisément ce qui mettra plus loin la validation croisée aléatoire en défaut.
Échauffement : ajuster la vitesse séculaire. Avant toute validation croisée, ajustons une droite sur toute la série, comme le ferait un géodésien pour estimer la vitesse de plaque.
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 pente ajustée s’écarte de la valeur vraie de 12 mm/an. Le saut cosismique et la décroissance postsismique tirent la droite vers le haut, et le bruit coloré ajoute une errance de longue période qu’une droite absorbe en partie. Un modèle peut être simple, bien ajusté, et néanmoins biaisé quand la physique contenue dans les données est plus riche que le modèle.
2. Une tâche de prévision supervisée¶
Les questions de validation croisée deviennent tranchantes quand le modèle fait quelque chose d’opérationnel. En voici une : étant donné l’histoire récente de la station, où sera-t-elle le mois prochain ? Pour chaque jour d’ancrage, nous construisons des caractéristiques à partir des 60 derniers jours et prédisons le déplacement moyen sur les 30 jours suivants.
Caractéristiques par jour d’ancrage :
- moyenne des 60 derniers jours
- moyenne des 10 derniers jours
- dernière valeur observée
- pente linéaire sur les 60 derniers jours (
np.polyfit) - sinus et cosinus du jour de l’année (pour encoder la saison)
Cible : déplacement moyen sur les 30 jours suivants.
Nous faisons glisser le jour d’ancrage avec un pas de 5 jours, ce qui donne environ 700 échantillons. Notez ce que fait cette construction : des jours d’ancrage voisins partagent presque toutes les données de leur fenêtre. L’échantillon i et l’échantillon i+1 se recouvrent sur 55 des 60 jours d’histoire et sur 25 des 30 jours de cible. Les échantillons sont fortement corrélés par construction, exactement comme ils le seraient pour n’importe quelle série temporelle géophysique fenêtrée.
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
Le modèle est le même pour toutes les expériences ci-dessous — une forêt aléatoire —, si bien que toute différence de score vient du découpage, pas du modèle.
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor(n_estimators=200, random_state=42)3. Le mensonge optimiste : la validation croisée aléatoire¶
D’abord, la recette standard qu’enseignent la plupart des tutoriels : mélanger les échantillons en 5 plis et noter avec cross_val_score. Nous essayons deux schémas mélangés, KFold(shuffle=True) et ShuffleSplit, et rapportons le R² et l’erreur absolue moyenne.
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
Les scores ont l’air excellents : R² proche de 1, erreurs d’un ou deux millimètres. Si c’était un article, le modèle serait déclaré réussi. Il ne l’est pas. La section suivante montre ce que note le même modèle quand le découpage respecte le temps.
4. Des découpages honnêtes pour les séries temporelles¶
TimeSeriesSplit entraîne sur le passé et valide sur le futur, ce qui correspond à l’usage réel d’un modèle de prévision. Comme cas intermédiaire, nous essayons aussi KFold sans mélange, qui valide sur des blocs de temps contigus : pas de quasi-doublons de part et d’autre de la frontière, sauf aux bords des blocs.
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
Même modèle, mêmes données, et les schémas honnêtes rapportent un skill (score de compétence) bien moindre. Pour comprendre pourquoi, traçons quels échantillons tombent en entraînement et lesquels en validation pour chaque schéma. L’indice d’échantillon court de gauche à droite, ce qui ici est aussi le temps.
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()
Dans les deux panneaux mélangés, les échantillons de validation en orange sont saupoudrés parmi les échantillons d’entraînement en bleu. Chaque échantillon de validation a des voisins immédiats dans l’ensemble d’entraînement, et ces voisins partagent 55 de ses 60 jours d’histoire. La forêt n’a rien besoin de prévoir : elle interpole ses voisins, et le score mesure de la mémorisation. TimeSeriesSplit (panneau du bas) valide toujours sur des données postérieures à tout ce sur quoi il a été entraîné. Le modèle doit extrapoler vers un futur qu’il n’a jamais vu, et le score reflète le skill qu’il aurait en exploitation.
Les blocs contigus (KFold sans mélange) se placent entre les deux : la fuite de données ne se produit qu’au voisinage des bords de blocs, si bien que les scores redescendent presque jusqu’au chiffre honnête. Les découpages par blocs sont la bonne idée dès qu’il n’y a pas d’axe temporel unique — par exemple des blocs spatiaux de stations.
5. Tout score exige un modèle de référence et une barre d’erreur¶
Les scores honnêtes de la section 4 ont l’air mauvais. À quel point ? Deux habitudes remettent n’importe quel score en contexte : le comparer à un modèle de référence (baseline) qui n’exige aucun apprentissage, et lui associer une barre d’erreur.
Le modèle de référence naturel pour une prévision est la persistance : le futur ressemble au présent. Ici, la prévision par persistance de la moyenne des 30 jours suivants est simplement la dernière valeur observée, une caractéristique que nous avons déjà calculée (last_val). Elle ne coûte rien et encode la propriété la plus forte de la série : elle change lentement.
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 comparaison est humiliante. La persistance prévoit la moyenne à 30 jours à environ 2 mm près et obtient un R² proche de 0,9 sur chaque pli ; la forêt est un ordre de grandeur moins bonne. La raison est l’extrapolation : la station continue de se déplacer à 12 mm/an, si bien que chaque cible de validation se situe au-dessus de la plage des cibles sur lesquelles la forêt s’est entraînée, et une forêt ne peut pas prédire hors de la plage de ses étiquettes d’entraînement. La persistance, elle, suit la tendance gratuitement.
Notez les valeurs de R² négatives du modèle. Le R² compare un modèle à la constante qui prédit la moyenne des cibles de validation, donc R² < 0 signifie que le modèle est pire que prédire la moyenne — aucun skill. Chaque fois qu’un article rapporte un skill de prévision sans modèle de référence de persistance (ou de climatologie), demandez-vous ce qu’aurait obtenu ce modèle de référence.
La barre d’erreur. Un chiffre unique comme « MAE = 16 mm » masque l’incertitude sur l’estimation elle-même. Le bootstrap lui associe un intervalle : rééchantillonner les erreurs de validation avec remise, recalculer la MAE de chaque rééchantillon, et lire un intervalle de confiance dans la dispersion. Pour des erreurs indépendantes, rééchantillonner des erreurs isolées (le bootstrap simple) suffit. Nos erreurs sont autocorrélées — des jours d’ancrage voisins partagent l’essentiel de leur fenêtre —, si bien que le rééchantillonnage par erreurs isolées fait comme s’il y avait plus d’information indépendante qu’il n’y en a, et l’intervalle ressort trop étroit. Le bootstrap par blocs mobiles rééchantillonne des blocs contigus (ici 12 échantillons, soit 60 jours au pas de 5 jours), de sorte que chaque rééchantillon conserve la corrélation locale.
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
L’intervalle par blocs est environ deux fois plus large que l’intervalle simple. Mêmes données, même statistique ; la seule différence est l’hypothèse d’indépendance. Des erreurs corrélées portent moins d’échantillons effectifs que ne le suggère leur nombre, et l’intervalle étroit du bootstrap simple est une forme de plus de l’optimisme dont parle toute cette leçon. Rapportez l’intervalle qui correspond à la structure de corrélation des erreurs, et rapportez-le à côté de l’estimation ponctuelle.
6. Validation croisée « laisser-un-dehors » (LOOCV)¶
La LOOCV découpe les données en ensembles d’entraînement et de validation n fois, où n est le nombre de points de données. À chaque tour, l’ensemble d’entraînement est constitué de tous les échantillons sauf un, et l’ensemble de validation de ce seul échantillon tenu à l’écart.
Avantages : faible biais vis-à-vis des données d’entraînement, et résultat déterministe, puisqu’il n’y a pas de découpage aléatoire à répéter. Inconvénient : elle coûte n ajustements de modèle, ce qui est cher pour tout ce qui n’est pas un petit jeu de données et un modèle bon marché.
Nous faisons la démonstration sur le problème simple de régression de la vitesse, sous-échantillonné à un jour sur 20 pour que les n ajustements restent rapides.
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 vitesse bouge à peine lorsqu’un point est retiré, si bien que la dispersion sur les n ajustements est minuscule. Deux mises en garde. D’abord, la LOOCV vaut rarement son coût : n ajustements pour une estimation d’erreur que la validation croisée à k plis approche avec 5 ou 10. Ensuite, la LOOCV ne corrige pas la corrélation. Les voisins immédiats de chaque point tenu à l’écart sont toujours dans l’ensemble d’entraînement, donc pour une série autocorrélée la LOOCV est proche du découpage le plus optimiste possible : c’est le problème des plis mélangés poussé à sa limite.
7. Recherche d’hyperparamètres sous validation croisée honnête¶
Les modèles ont des paramètres appris des données (coupures des arbres, poids de régression) et des hyperparamètres fixés avant l’entraînement (profondeur des arbres, taille des feuilles). Le réglage des hyperparamètres cherche les réglages qui obtiennent le meilleur score en validation croisée ; c’est une pratique standard. Les approches usuelles :
- Réglage manuel : ajuster à la main, guidé par la connaissance du problème. Bon pour construire l’intuition ; pas systématique.
- Recherche exhaustive sur grille : évaluer chaque combinaison d’une grille prédéfinie. Complète mais coûteuse à mesure que la grille grandit. Dans scikit-learn :
GridSearchCV. - Recherche aléatoire : tirer des combinaisons dans des distributions prédéfinies, pour un budget fixé d’itérations. Couvre de larges espaces à moindre coût qu’une grille. Dans scikit-learn :
RandomizedSearchCV. - Optimisation bayésienne : modéliser le score comme une fonction des hyperparamètres et dépenser les évaluations là où une amélioration semble probable (p. ex.
optuna,scikit-optimize).
Le choix du schéma de validation croisée à l’intérieur de la recherche compte autant que la recherche elle-même. Régler sur des plis mélangés sélectionne le modèle qui mémorise le mieux ses voisins. Nous réglons sur TimeSeriesSplit, si bien que le vainqueur est le meilleur prévisionniste.
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
Pour la recherche aléatoire, passez les objets de distribution figés (frozen) de scipy.stats. Une version antérieure de ce carnet pré-échantillonnait avec randint.rvs(size=10), ce qui réduit la recherche aléatoire à une liste fixe de dix valeurs ; passer la distribution figée permet à chaque itération de tirer de nouvelles valeurs.
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
Les deux recherches arrivent à une fraction de millimètre l’une de l’autre : avec seulement deux hyperparamètres, 15 tirages aléatoires sondent l’espace à peu près aussi bien qu’une grille à 12 points, et la recherche aléatoire passe mieux à l’échelle quand l’espace grandit.
8. Espace et groupes : le même mensonge sur une carte¶
Le temps n’est pas le seul axe le long duquel les données géoscientifiques sont corrélées. Les campagnes de terrain produisent des données groupées : plusieurs sondages par site, plusieurs sites le long de chaque route ou de chaque bassin, et une cible qui repose sur un champ régional lisse. Deux observations d’un même site sont de quasi-copies ; deux sites d’un même groupe spatial sont de proches cousins.
mlgeo_synth.multisite_table construit cette géométrie avec une vérité terrain : 6 groupes de 5 sites chacun, 10 observations répétées par site (300 lignes). La cible est un signal linéaire dans les colonnes feat_* — la part transportable du skill — plus un champ régional lisse (longueur de corrélation de 25 km) évalué à chaque site, plus du bruit. Le champ est renvoyé sous forme d’objet appelable dans truth, ce qui nous permet de dessiner la carte que le modèle mémorise en douce.
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()
La carte est le mécanisme. Le champ varie doucement sur ~25 km, et les sites se tiennent en groupes serrés de quelques kilomètres de large, si bien que tous les sites d’un groupe partagent presque la même valeur de champ — et les 10 observations d’un site la partagent exactement. Un découpage mélangé livre au modèle ces quasi-copies.
Voici maintenant l’échelle des découpages. Même forêt, mêmes données, trois découpages qui tiennent à l’écart une part progressivement plus grande de la structure partagée : KFold mélangé, GroupKFold sur site_id (sites entiers tenus à l’écart), et GroupKFold sur cluster_id avec 6 découpages — laisser-un-groupe-dehors (leave-one-cluster-out), la version discrète d’un découpage spatial par blocs. Nous parcourons cette échelle deux fois : une fois sur les seules colonnes feat_*, et une fois en ajoutant les coordonnées x_km, y_km comme caractéristiques.
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)Lisez l’échelle de haut en bas ; chaque barreau pose une question plus difficile, et plus honnête.
- KFold (mélangé), R² ≈ 0,7. Chaque observation tenue à l’écart a jusqu’à 9 sœurs de son propre site dans l’entraînement. La forêt reconnaît le site à ses caractéristiques et se rappelle la valeur de champ de ce site. Le score répond à : « avec quelle qualité puis-je prédire un autre sondage sur un site que j’ai déjà prospecté ? »
- GroupKFold par site, R² ≈ 0,35. Des sites entiers sont tenus à l’écart, donc la mémorisation du site disparaît — mais les sites tenus à l’écart ont encore, dans l’entraînement, des voisins du même groupe qui partagent le champ. Le score répond à : « un nouveau site à l’intérieur d’un groupe prospecté ? »
- Laisser-un-groupe-dehors, R² < 0. Des groupes entiers sont tenus à l’écart et le modèle affronte une région qu’il n’a jamais vue. La valeur du champ y est inconnaissable à partir des données d’entraînement, et les décalages de champ mémorisés par la forêt l’induisent activement en erreur — pire que prédire la moyenne (section 5). Le score répond à la question que pose réellement une carte régionale d’aléa : « une nouvelle région ? »
Ajouter x_km, y_km accentue le contraste. Sous validation croisée mélangée, la forêt interpole désormais la carte presque parfaitement (R² ≈ 0,95) — un skill authentiquement utile à l’intérieur de la région prospectée. Sous laisser-un-groupe-dehors, les mêmes coordonnées aggravent les choses : elles pointent vers des parties de la carte que les données d’entraînement n’ont jamais contraintes. Même modèle, mêmes caractéristiques, verdicts opposés — parce qu’interpoler entre des sites prospectés et extrapoler vers une nouvelle région sont deux affirmations différentes, et que le découpage décide laquelle des deux le score étaie. Dans des contextes continus sans groupes naturels, la même idée devient un découpage spatial avec zone tampon : exclure de l’entraînement tout ce qui se trouve à moins d’une longueur de corrélation des sites de validation.
Le carnet 4.5 (section 4.6) franchit cette ligne délibérément : un ensemble profond entraîné sur une plage de compositions bornée y est évalué hors plage — y compris sur la seule erreur que son désaccord ne parvient pas à signaler.
9. Petit, groupé dans l’espace, déséquilibré : quand la validation croisée par groupes se casse elle-même¶
En pratique géotechnique et géologique, les données groupées sont généralement aussi petites et déséquilibrées : quelques centaines de cas historiques, une classe positive rare (liquéfaction observée, glissement de terrain survenu, minéralisation présente), et des positifs qui se regroupent dans l’espace parce que le champ qui les pilote le fait aussi. Cette combinaison est la norme, pas le cas particulier, partout où les étiquettes proviennent de cas historiques plutôt que de levés. multisite_table(binary=True) la reproduit : la même géométrie de sites à 300 lignes, avec la cible latente seuillée de sorte que 12 % des lignes soient positives.
Grouper par site reste obligatoire — mais peut désormais échouer selon ses propres termes. GroupKFold distribue des sites entiers sans regarder les étiquettes, et comme les positifs sont regroupés spatialement, un pli peut se retrouver sans aucun positif. L’AUC ROC est indéfinie sur un pli qui ne contient qu’une classe.
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]
Le nan n’est pas un bogue de scikit-learn ; ce sont les données qui signalent que ce découpage par groupes a laissé un pli avec zéro positif, si bien qu’il n’y a pas de courbe ROC à calculer — et moyenner les plis restants change discrètement ce que la moyenne estime. La réparation, c’est StratifiedGroupKFold : garder chaque site intact et équilibrer la fraction de positifs entre les plis.
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
Chaque pli a maintenant des positifs et chaque pli renvoie un nombre. Conservez la dispersion d’un pli à l’autre dans le rapport : avec 36 positifs répartis en cinq, l’AUC de chaque pli repose sur une poignée d’événements, et la dispersion (ou un intervalle bootstrap, section 5) est aussi informative que la moyenne.
10. Grouper par événement : répliques et tirs de carrière¶
En sismologie, le groupe est en général l’événement. mlgeo_synth.event_station_table construit un jeu de données de mouvement du sol jouet : 60 séismes en 4 groupes (compacts dans l’espace et dans le temps, comme des séquences choc principal–répliques), chacun enregistré par les mêmes 15 stations — 900 enregistrements. Le logarithme de l’accélération maximale du sol suit une relation d’atténuation jouet en magnitude et en distance, plus un terme par événement partagé par les 15 enregistrements d’un événement et un terme par station. Le terme d’événement est la fuite : un découpage aléatoire met 12 enregistrements d’un événement en entraînement et 3 en validation, et le modèle est crédité pour s’être rappelé le décalage de cet événement.
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
Grouper par événement coûte environ un quart du skill apparent, et c’est le chiffre groupé qui prédit la performance sur le prochain séisme. La version classification de ce piège est fréquente : des centaines de répliques d’un même choc principal, ou des tirs répétés d’une même carrière, sont de quasi-copies les unes des autres, et un découpage aléatoire permet à un classifieur d’obtenir un score élevé en reconnaissant la source plutôt que le type de source. C’est exactement la fuite que l’encadré d’avertissement de 3.5 déclare et tolère — les tableaux de caractéristiques du classement (leaderboard) ne portent aucun identifiant d’événement sur lequel grouper. Ici, les métadonnées existent, donc le coût de cette concession est mesurable : c’est l’écart entre les deux lignes ci-dessus.
11. Choisir le découpage¶
Chaque schéma de cette leçon répond à la même question, posée à une échelle différente de structure partagée. Avant de faire confiance à un score, demandez-vous : quelle structure mes données partagent-elles, que le découpage doit respecter — le temps, le site, l’événement ou l’espace ? Des observations partageant une fenêtre de temps appellent TimeSeriesSplit ou des blocs contigus ; partageant un site ou un instrument, GroupKFold par site ; partageant un événement source, GroupKFold par événement ; partageant une région, laisser-un-groupe-dehors ou un découpage spatial avec zone tampon. L’audit de fuite de données de 2.13 pose cette question avant l’entraînement ; le choix du découpage est le même audit appliqué à l’évaluation. Et quand les métadonnées nécessaires au groupement manquent, dites-le et énoncez ce que le score peut et ne peut pas signifier, comme le fait l’encadré du classement de 3.5.