Prévoir, c’est prédire les valeurs futures d’une série temporelle à partir de son passé. Les géosciences fonctionnent à la prévision : CO2 atmosphérique, débit des rivières, écoulement des glaciers, niveaux des nappes, météorologie de l’espace. Ce carnet est un concours de prévision. Nous prenons les familles de modèles canoniques rencontrées dans ce livre, nous les pointons toutes vers les deux mêmes jeux de données réels, nous les évaluons avec les mêmes métriques sur le même découpage temporel, et nous laissons un tableau comparatif trancher.
Les concurrents, par complexité croissante :
- Les références (baselines) : persistance naïve et naïf saisonnier. Toujours en premier.
- SARIMA : un modèle statistique classique (
statsmodels). - Le gradient boosting sur caractéristiques de retard (
lightgbm) : l’apprentissage automatique classique du chapitre 3, adapté au temps. - Un petit LSTM (
torch) : l’architecture récurrente du carnet 4.4. - Un petit encodeur transformer (
torch) : le motif d’attention du carnet 4.4.
Les données :
- CO2 mensuel de Mauna Loa (NOAA Global Monitoring Laboratory) : longue, régulière, fortement saisonnière, avec une tendance lisse. Le cas facile, et notre référence principale.
- Vitesse de surface du glacier Jakobshavn Isbræ (Groenland) : courte, échantillonnée irrégulièrement, bruitée. Le cas difficile.
Objectifs d’apprentissage
- Découper des séries temporelles dans le temps et expliquer pourquoi les découpages aléatoires laissent fuiter de l’information future.
- Construire des caractéristiques de retard qui n’utilisent que le passé.
- Évaluer des prévisions avec MSE, RMSE, MAE, MAPE et MASE, calculées dans le code.
- Transformer un prévisionniste ponctuel en prévisionniste probabiliste avec la perte quantile (pinball loss), puis le vérifier : couverture empirique de l’intervalle à 90 % et CRPS.
- Mesurer comment l’erreur de prévision croît avec l’échéance et localiser l’horizon de skill (skill horizon, horizon du score de compétence).
- Quantifier de combien un classement de modèles bouge d’une graine aléatoire à l’autre, et noter les événements rares séparément du gros de la série.
- Juger quand la complexité d’un modèle est rentable et quand une référence l’emporte.
Le carnet se termine par un classement (leaderboard) de la classe en deux pistes : une prévision publique de CO2 qui n’est qu’un diagnostic (la section 7 explique pourquoi elle est trivialement contournable, et pourquoi nous l’écrivons noir sur blanc) et une série synthétique cachée qui porte tout le poids de la note.
import os
import warnings
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pooch
import torch
import torch.nn as nn
import mlgeo_synth
np.random.seed(0)
torch.manual_seed(0)
device = torch.device("cuda" if torch.cuda.is_available()
else "mps" if torch.backends.mps.is_available()
else "cpu")
print("device:", device)
# One fixed color per model, reused in every figure (Okabe-Ito palette, colorblind-safe).
COLORS = {
"observations": "#000000",
"naive": "#999999",
"seasonal naive": "#E69F00",
"SARIMA": "#0072B2",
"ARIMA": "#0072B2",
"LightGBM": "#009E73",
"LSTM": "#D55E00",
"Transformer": "#CC79A7",
}
plt.rcParams.update({"figure.dpi": 100, "axes.grid": True, "grid.alpha": 0.3})device: mps
1. Les données¶
1.1 CO2 mensuel de Mauna Loa¶
L’enregistrement de Mauna Loa est la plus longue mesure directe continue du CO2 atmosphérique, entamée par Charles Keeling en 1958. La NOAA GML en distribue les moyennes mensuelles dans un fichier texte : environ 42 lignes de commentaires #, une ligne d’en-tête, puis des colonnes séparées par des espaces (année, mois, date décimale, moyenne mensuelle en ppm, valeur désaisonnalisée, nombre de jours, écart-type, incertitude). Nous le récupérons avec pooch, qui met le fichier en cache localement.
co2_url = "https://raw.githubusercontent.com/UW-MLGEO/MLGeo-dataset/main/data/data_co2.csv"
co2_file = pooch.retrieve(co2_url, known_hash=None, fname="data_co2.csv")
co2 = pd.read_csv(co2_file, comment="#", sep=r"\s+", header=0)
co2.columns = ["year", "month", "decimal_date", "average",
"deseasonalized", "ndays", "stdev", "unc"]
co2.index = pd.to_datetime(co2["year"].astype(str) + "-"
+ co2["month"].astype(str).str.zfill(2))
y_full = co2["average"].astype("float64")
print(y_full.shape, "monthly values,", y_full.index[0].date(), "to", y_full.index[-1].date())(798,) monthly values, 1958-03-01 to 2024-08-01
1.2 Les données réservées du classement : bas les mains sur les 12 derniers mois¶
Avant toute modélisation, nous retirons les 12 derniers mois de l’enregistrement : de 2023-09 à 2024-08. Ces mois constituent l’ensemble réservé du classement de la classe (section 7) — un ensemble de test que nous noterons exactement une fois. Rien dans les sections de modélisation n’y touche : ni figure, ni ajustement, ni métrique. Nous n’en gardons que les dates.
holdout_months = y_full.index[-12:]
y = y_full.iloc[:-12].copy()
del y_full # the holdout values are gone from this notebook
assert holdout_months[0] == pd.Timestamp("2023-09-01")
assert holdout_months[-1] == pd.Timestamp("2024-08-01")
print("holdout months:", holdout_months.strftime("%Y-%m").tolist())
print("working series:", y.index[0].date(), "to", y.index[-1].date(), f"({len(y)} months)")holdout months: ['2023-09', '2023-10', '2023-11', '2023-12', '2024-01', '2024-02', '2024-03', '2024-04', '2024-05', '2024-06', '2024-07', '2024-08']
working series: 1958-03-01 to 2023-08-01 (786 months)
fig, axes = plt.subplots(1, 2, figsize=(10, 3.2))
axes[0].plot(y.index, y.values, color=COLORS["observations"], lw=1)
axes[0].set_title("Mauna Loa CO$_2$, full record")
axes[1].plot(y.loc["2013":].index, y.loc["2013":].values,
color=COLORS["observations"], lw=1.2)
axes[1].set_title("Last decade")
for ax in axes:
ax.set_xlabel("year")
ax.set_ylabel("CO$_2$ (ppm)")
fig.tight_layout()
plt.show()
Une tendance lisse et accélérée, plus un cycle annuel en dents de scie (la végétation de l’hémisphère nord pompe le CO2 chaque été). Une structure hautement prévisible. Tout modèle incapable de battre « répéter la forme de l’an dernier » devrait ici faire honte à son auteur.
1.3 Vitesse de surface du Jakobshavn Isbræ¶
Le Jakobshavn Isbræ, dans l’ouest du Groenland, est l’un des glaciers les plus rapides de la planète. Le fichier contient la vitesse de surface dérivée de données satellitaires (m/an) pour une imagette de 10 x 10 pixels près du front, échantillonnée irrégulièrement (tous les 6 à 12 jours) depuis 2015. Les valeurs manquantes sont codées -1.0 ; nous les masquons en NaN et moyennons les pixels valides à chaque date en une seule série de vitesse moyenne.
ice_url = "https://raw.githubusercontent.com/UW-MLGEO/MLGeo-dataset/main/data/data_ice_jakobshavn.csv"
ice_file = pooch.retrieve(ice_url, known_hash=None, fname="data_ice_jakobshavn.csv")
ice_raw = pd.read_csv(ice_file, parse_dates=["Date"]).set_index("Date")
ice_raw = ice_raw.where(ice_raw != -1.0) # -1.0 means missing
ice = ice_raw.mean(axis=1).dropna() # patch mean over valid pixels
print(ice.shape, "samples,", ice.index[0].date(), "to", ice.index[-1].date())
print("median sampling interval:", ice.index.to_series().diff().median())
fig, ax = plt.subplots(figsize=(9, 3))
ax.plot(ice.index, ice.values, color=COLORS["observations"], lw=1)
ax.set_xlabel("year")
ax.set_ylabel("surface speed (m/yr)")
ax.set_title("Jakobshavn Isbrae, patch-mean surface speed")
fig.tight_layout()
plt.show()(430,) samples, 2015-01-01 to 2024-04-13
median sampling interval: 6 days 00:00:00
Les accélérations estivales se superposent à des oscillations pluriannuelles liées à la position du front et au forçage océanique. Comparée au CO2, cette série est courte (environ 430 échantillons sur neuf ans), échantillonnée irrégulièrement, et sa variabilité n’est pas dominée par un cycle répétitif net.
2. Protocole d’évaluation¶
2.1 Découpages temporels, identiques pour tous les modèles¶
Un modèle de prévision n’est équitablement testé que sur des données postérieures à tout ce sur quoi il a été entraîné. Un découpage entraînement/test aléatoire des lignes d’une série temporelle laisse le modèle interpoler entre les voisins des points de test, ce qui gonfle les scores et ne prévoit rien. Tous les modèles ci-dessous emploient exactement ces découpages :
- CO2 : TEST = les 48 derniers mois de la série de travail, de 2019-09 à 2023-08. TRAIN = tout ce qui précède, de 1958-03 à 2019-08.
- Glace : TEST = les 20 % de dates les plus récentes. TRAIN = les 80 % les plus anciennes.
Tous les modèles prévoient la fenêtre de test entière d’un seul coup (un horizon de 48 mois pour le CO2), sans jamais jeter un œil aux valeurs de test en chemin.
H_CO2 = 48
train_co2, test_co2 = y.iloc[:-H_CO2], y.iloc[-H_CO2:]
assert test_co2.index[0] == pd.Timestamp("2019-09-01")
assert test_co2.index[-1] == pd.Timestamp("2023-08-01")
n_test_ice = int(0.2 * len(ice))
train_ice, test_ice = ice.iloc[:-n_test_ice], ice.iloc[-n_test_ice:]
print(f"CO2 : train {train_co2.index[0].date()} .. {train_co2.index[-1].date()} "
f"({len(train_co2)}), test {test_co2.index[0].date()} .. {test_co2.index[-1].date()} ({len(test_co2)})")
print(f"ice : train {train_ice.index[0].date()} .. {train_ice.index[-1].date()} "
f"({len(train_ice)}), test {test_ice.index[0].date()} .. {test_ice.index[-1].date()} ({len(test_ice)})")CO2 : train 1958-03-01 .. 2019-08-01 (738), test 2019-09-01 .. 2023-08-01 (48)
ice : train 2015-01-01 .. 2021-09-20 (344), test 2021-09-26 .. 2024-04-13 (86)
2.2 Les métriques¶
Avec les erreurs de prévision sur les points de test :
- MSE . Pénalise quadratiquement les grandes erreurs ; ses unités sont au carré, donc malcommodes à lire.
- RMSE . De retour dans les unités des données ; toujours dominée par les plus grandes erreurs.
- MAE . Unités des données, toutes les erreurs traitées linéairement. La plus facile à interpréter.
- MAPE . Un pourcentage sans unité, mais qui explose près des valeurs nulles et qui est asymétrique (sur-prévisions et sous-prévisions sont pénalisées différemment). Convenable pour le CO2 (~420 ppm) et la vitesse de la glace, inutilisable pour des séries qui traversent zéro.
- MASE (mean absolute scaled error, erreur absolue moyenne mise à l’échelle) : la MAE divisée par la MAE de la prévision naïve à un pas () sur les données d’entraînement. Une MASE signifie que le modèle bat la persistance naïve à son propre jeu à un pas ; une MASE signifie qu’une règle « répéter la dernière valeur » faisait mieux. Elle est sans unité, définie en zéro et comparable d’une série à l’autre, ce qui explique que le classement de la classe l’emploie.
Une seule fonction utilitaire les calcule toutes les cinq. Tous les nombres des tableaux comparatifs sortent de cette fonction. Ne saisissez jamais une valeur de métrique à la main dans un tableau ; l’édition 2024 de ce carnet l’a fait, et le tableau était faux.
def forecast_metrics(y_true, y_pred, y_train):
"""All metrics for one forecast. MASE is scaled by the training
MAE of the naive one-step forecast (y_hat_t = y_{t-1})."""
y_true = np.asarray(y_true, dtype=float)
y_pred = np.asarray(y_pred, dtype=float)
e = y_true - y_pred
mse = np.mean(e**2)
mae = np.mean(np.abs(e))
naive_train_mae = np.mean(np.abs(np.diff(np.asarray(y_train, dtype=float))))
return {"MSE": mse,
"RMSE": np.sqrt(mse),
"MAE": mae,
"MAPE (%)": 100 * np.mean(np.abs(e / y_true)),
"MASE": mae / naive_train_mae}
results_co2, preds_co2 = {}, {} # metrics and predictions, filled model by model
results_ice, preds_ice = {}, {}
def register(results, preds, name, y_pred, test, train):
preds[name] = np.asarray(y_pred, dtype=float)
results[name] = forecast_metrics(test.values, y_pred, train.values)
print(name, "->", {k: round(v, 3) for k, v in results[name].items()})3. Le concours sur le CO2¶
3.1 Les références d’abord, toujours¶
- Persistance naïve : répéter la dernière valeur d’entraînement sur tout l’horizon de 48 mois.
- Naïf saisonnier : répéter les 12 derniers mois d’entraînement, dupliqués quatre fois. Cela recopie la forme saisonnière de l’an dernier mais gèle la tendance.
Si un modèle ne peut pas battre ces références, sa complexité n’a rien acheté.
naive = np.full(H_CO2, train_co2.iloc[-1])
seasonal = np.tile(train_co2.values[-12:], H_CO2 // 12)
register(results_co2, preds_co2, "naive", naive, test_co2, train_co2)
register(results_co2, preds_co2, "seasonal naive", seasonal, test_co2, train_co2)naive -> {'MSE': np.float64(56.38), 'RMSE': np.float64(7.509), 'MAE': np.float64(6.679), 'MAPE (%)': np.float64(1.596), 'MASE': np.float64(6.09)}
seasonal naive -> {'MSE': np.float64(41.971), 'RMSE': np.float64(6.478), 'MAE': np.float64(5.961), 'MAPE (%)': np.float64(1.427), 'MASE': np.float64(5.436)}
3.2 SARIMA¶
SARIMA (ARIMA saisonnier) modélise la série différenciée avec des termes autorégressifs et de moyenne mobile, doublés de leurs équivalents saisonniers au retard 12. Nous employons l’ordre canonique des manuels pour Mauna Loa, : une différenciation simple pour la tendance, une différenciation saisonnière pour le cycle annuel. Nous plafonnons le nombre d’itérations de l’optimiseur pour que l’ajustement reste rapide ; c’est une entrée de concours, pas un exercice de réglage.
from statsmodels.tsa.statespace.sarimax import SARIMAX
with warnings.catch_warnings():
warnings.simplefilter("ignore")
sarima = SARIMAX(train_co2, order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12)).fit(disp=False, maxiter=50)
sarima_pred = sarima.forecast(H_CO2).values
register(results_co2, preds_co2, "SARIMA", sarima_pred, test_co2, train_co2)SARIMA -> {'MSE': np.float64(0.129), 'RMSE': np.float64(0.359), 'MAE': np.float64(0.281), 'MAPE (%)': np.float64(0.067), 'MASE': np.float64(0.256)}
3.3 Gradient boosting sur caractéristiques de retard¶
Les ensembles d’arbres ne savent rien du temps. Nous le leur donnons explicitement, sous forme d’un tableau où chaque ligne prédit un mois à partir de caractéristiques construites uniquement à partir des mois antérieurs :
- les retards 1, 2, 3, 6, 12, 24 (la valeur mois avant la cible),
- la moyenne glissante des 12 mois précédant la cible (c’est-à-dire décalée d’un pas, de sorte que la cible elle-même soit exclue),
- le mois calendaire (connu à l’avance, donc de bonne guerre).
Une subtilité : un arbre de régression prédit des moyennes de cibles d’entraînement, il ne peut donc jamais produire une valeur hors de la plage vue à l’entraînement. Or le CO2 continue de monter au-delà de toute valeur d’entraînement : un arbre qui prédit le niveau sature près du maximum d’entraînement et la prévision devient plate. Le remède consiste à prédire la variation mensuelle et à l’ajouter à la valeur précédente ; les variations mensuelles sont à peu près stationnaires et se situent bien à l’intérieur de la plage d’entraînement. Les caractéristiques restent les mêmes retards, uniquement passés.
Le modèle prédit un pas en avant. Pour couvrir la fenêtre de test de 48 mois, nous prévoyons récursivement : prédire la variation suivante, l’ajouter à la dernière valeur, ajouter cette prédiction à l’historique, reconstruire les caractéristiques, et ainsi de suite. Les erreurs se composent, ce qui est honnête ; le modèle ne reçoit aucune observation réelle de la période de test.
from lightgbm import LGBMRegressor
def lag_feature_row(history, date, lags, roll):
"""Features for predicting the value at `date`, given all values before it."""
row = {f"lag_{k}": history[-k] for k in lags}
row[f"rollmean_{roll}"] = np.mean(history[-roll:])
row["month"] = date.month
return row
def build_training_table(series, lags, roll):
"""One row per predictable month; the target is the one-step change."""
vals = series.to_numpy()
start = max(max(lags), roll) # first row with a full feature set
rows = [lag_feature_row(vals[:i], series.index[i], lags, roll)
for i in range(start, len(vals))]
target = vals[start:] - vals[start - 1:-1] # y_t - y_{t-1}
return pd.DataFrame(rows), target
def lgbm_recursive_forecast(train_series, test_index, lags, roll, random_state=0):
X_tr, y_tr = build_training_table(train_series, lags, roll)
model = LGBMRegressor(n_estimators=300, learning_rate=0.05,
num_leaves=15, random_state=random_state, verbose=-1)
model.fit(X_tr, y_tr)
history = list(train_series.to_numpy())
preds = []
for date in test_index: # recursive multi-step
row = pd.DataFrame([lag_feature_row(np.asarray(history), date, lags, roll)])
step = float(model.predict(row)[0]) # predicted change
level = history[-1] + step
preds.append(level)
history.append(level) # feed the prediction back in
return np.array(preds), modellgbm_pred_co2, _ = lgbm_recursive_forecast(train_co2, test_co2.index,
lags=[1, 2, 3, 6, 12, 24], roll=12)
register(results_co2, preds_co2, "LightGBM", lgbm_pred_co2, test_co2, train_co2)LightGBM -> {'MSE': np.float64(1.643), 'RMSE': np.float64(1.282), 'MAE': np.float64(1.03), 'MAPE (%)': np.float64(0.246), 'MASE': np.float64(0.939)}
3.4 Un petit LSTM¶
Les modèles profonds prévoient directement : une seule passe avant (forward pass) transforme une fenêtre de contexte de 36 mois en l’ensemble des 48 valeurs de l’horizon d’un coup (pas de récursion). Deux points pratiques :
- Ancrage. Le CO2 croît indéfiniment, si bien que les fenêtres de test se situent à des niveaux que le réseau n’a jamais vus à l’entraînement. Nous soustrayons la dernière valeur du contexte (l’« ancre ») de chaque fenêtre, nous entraînons sur les formes ancrées, et nous rajoutons l’ancre à la prédiction. Le réseau apprend la variation par rapport à maintenant, suffisamment stationnaire pour généraliser.
- Mise à l’échelle. Les fenêtres ancrées sont divisées par un unique écart-type global calculé sur les fenêtres d’entraînement, ce qui maintient les entrées en O(1).
Le modèle est petit et s’entraîne quelques époques ; augmentez les deux sur votre propre machine si vous le souhaitez.
def make_windows(values, context, horizon):
"""Anchored (context, horizon) windows: both parts minus the last context value."""
X, Y = [], []
for i in range(len(values) - context - horizon + 1):
c = values[i:i + context]
h = values[i + context:i + context + horizon]
X.append(c - c[-1])
Y.append(h - c[-1])
return np.array(X, dtype=np.float32), np.array(Y, dtype=np.float32)
def train_forecaster(model, X, Y, epochs=40, lr=1e-3, batch_size=32):
model.to(device)
ds = torch.utils.data.TensorDataset(torch.from_numpy(X).unsqueeze(-1),
torch.from_numpy(Y))
dl = torch.utils.data.DataLoader(ds, batch_size=batch_size, shuffle=True)
opt = torch.optim.Adam(model.parameters(), lr=lr)
loss_fn = nn.MSELoss()
losses = []
for _ in range(epochs):
model.train()
total = 0.0
for xb, yb in dl:
xb, yb = xb.to(device), yb.to(device)
opt.zero_grad()
loss = loss_fn(model(xb), yb)
loss.backward()
opt.step()
total += loss.item() * len(xb)
losses.append(total / len(ds))
return losses
def nn_forecast(model, context_values, scale):
"""Forecast one horizon from the last `context` observed values."""
anchor = context_values[-1]
x = torch.from_numpy(((context_values - anchor) / scale)
.astype(np.float32)).reshape(1, -1, 1).to(device)
model.eval()
with torch.no_grad():
out = model(x).cpu().numpy().ravel()
return out * scale + anchorCONTEXT_CO2 = 36
Xc, Yc = make_windows(train_co2.to_numpy(), CONTEXT_CO2, H_CO2)
scale_co2 = float(Xc.std())
Xc, Yc = Xc / scale_co2, Yc / scale_co2
print(f"{len(Xc)} training windows, scale = {scale_co2:.2f} ppm")
class LSTMForecaster(nn.Module):
def __init__(self, horizon, hidden=32):
super().__init__()
self.lstm = nn.LSTM(input_size=1, hidden_size=hidden, batch_first=True)
self.head = nn.Linear(hidden, horizon)
def forward(self, x):
_, (h, _) = self.lstm(x)
return self.head(h[-1])
torch.manual_seed(0)
lstm_co2 = LSTMForecaster(horizon=H_CO2)
lstm_losses = train_forecaster(lstm_co2, Xc, Yc, epochs=40)
lstm_pred_co2 = nn_forecast(lstm_co2, train_co2.to_numpy()[-CONTEXT_CO2:], scale_co2)
register(results_co2, preds_co2, "LSTM", lstm_pred_co2, test_co2, train_co2)655 training windows, scale = 3.35 ppm
LSTM -> {'MSE': np.float64(1.085), 'RMSE': np.float64(1.042), 'MAE': np.float64(0.886), 'MAPE (%)': np.float64(0.212), 'MASE': np.float64(0.808)}
3.5 Un petit encodeur transformer¶
Mêmes fenêtres, même ancrage, même tête directe multi-pas ; seul le modèle de séquence change. Nous réutilisons le motif du carnet 4.4 : projeter chaque scalaire dans un petit plongement, ajouter un encodage positionnel sinusoïdal (l’attention est en elle-même aveugle à l’ordre), enchaîner deux couches d’encodeur, moyenner dans le temps, et projeter vers les 48 valeurs de l’horizon.
class PositionalEncoding(nn.Module):
def __init__(self, d_model, max_len=200):
super().__init__()
pos = torch.arange(max_len).unsqueeze(1).float()
div = torch.exp(torch.arange(0, d_model, 2).float()
* (-np.log(10000.0) / d_model))
pe = torch.zeros(max_len, d_model)
pe[:, 0::2] = torch.sin(pos * div)
pe[:, 1::2] = torch.cos(pos * div)
self.register_buffer("pe", pe)
def forward(self, x): # x: (batch, seq, d_model)
return x + self.pe[: x.shape[1]]
class TransformerForecaster(nn.Module):
def __init__(self, horizon, d_model=32, nhead=4, num_layers=2):
super().__init__()
self.embed = nn.Linear(1, d_model)
self.pos = PositionalEncoding(d_model)
layer = nn.TransformerEncoderLayer(d_model=d_model, nhead=nhead,
dim_feedforward=64, dropout=0.1,
batch_first=True)
self.encoder = nn.TransformerEncoder(layer, num_layers=num_layers)
self.head = nn.Linear(d_model, horizon)
def forward(self, x):
z = self.encoder(self.pos(self.embed(x)))
return self.head(z.mean(dim=1)) # average over time steps
torch.manual_seed(0)
tfm_co2 = TransformerForecaster(horizon=H_CO2)
tfm_losses = train_forecaster(tfm_co2, Xc, Yc, epochs=40)
tfm_pred_co2 = nn_forecast(tfm_co2, train_co2.to_numpy()[-CONTEXT_CO2:], scale_co2)
register(results_co2, preds_co2, "Transformer", tfm_pred_co2, test_co2, train_co2)Transformer -> {'MSE': np.float64(0.214), 'RMSE': np.float64(0.463), 'MAE': np.float64(0.358), 'MAPE (%)': np.float64(0.086), 'MASE': np.float64(0.327)}
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(lstm_losses, color=COLORS["LSTM"], lw=1.5, label="LSTM")
ax.plot(tfm_losses, color=COLORS["Transformer"], lw=1.5, label="Transformer")
ax.set_xlabel("epoch")
ax.set_ylabel("training MSE (scaled units)")
ax.set_yscale("log")
ax.set_title("Training loss, CO$_2$ deep models")
ax.legend()
fig.tight_layout()
plt.show()
3.6 La comparaison, calculée honnêtement¶
Le tableau ci-dessous est construit programmatiquement à partir de results_co2, le dictionnaire que register() a rempli avec les prédictions réelles. Aucun nombre n’y a été saisi à la main.
co2_table = pd.DataFrame(results_co2).T.sort_values("MASE").round(3)
co2_tablefig, ax = plt.subplots(figsize=(9, 4))
recent = y.loc["2017":]
ax.plot(recent.index, recent.values, color=COLORS["observations"],
lw=1.5, label="observations")
for name, pred in preds_co2.items():
ax.plot(test_co2.index, pred, color=COLORS[name], lw=1.2, ls="--", label=name)
ax.axvline(test_co2.index[0], color="0.6", lw=0.8)
ax.text(test_co2.index[0], ax.get_ylim()[1], " test window starts",
va="top", fontsize=8, color="0.4")
ax.set_xlabel("year")
ax.set_ylabel("CO$_2$ (ppm)")
ax.set_title("CO$_2$ test window: all forecasts vs observations")
ax.legend(fontsize=8, ncol=2)
fig.tight_layout()
plt.show()
Ce qu’il faut chercher dans votre tableau et votre figure :
- La persistance naïve ignore à la fois la tendance et la saison ; sa ligne plate est le plancher que tout le monde devrait battre, et sa MASE est grande parce qu’un horizon de 48 mois est bien plus difficile que le jeu à un pas par lequel la MASE est mise à l’échelle.
- Le naïf saisonnier recopie la forme annuelle mais gèle la tendance : il dérive donc sous les observations d’environ la tendance multipliée par l’horizon. Il bat malgré tout plusieurs modèles entraînés. C’est la leçon récurrente des concours de prévision.
- SARIMA encode exactement les deux ingrédients dont cette série est faite, une tendance différenciée et un cycle de 12 mois, en une poignée de paramètres. Sur ce genre de série, il est très difficile à battre.
- LightGBM, le LSTM et le transformer doivent apprendre tendance et saison à partir d’exemples. Dans une exécution typique, tous trois atteignent une MASE inférieure à 1 et se placent entre la référence saisonnière et SARIMA ; la rétroaction récursive (LightGBM) et le peu de données d’entraînement (les deux modèles profonds) coûtent l’une et l’autre en exactitude sur un horizon de 4 ans. Le goulot d’étranglement n’est pas la capacité du modèle, c’est la structure du problème.
3.7 De la prévision ponctuelle aux intervalles : perte quantile, couverture et CRPS¶
Tous les modèles vus jusqu’ici répondent « quel sera le CO2 ? » par un seul nombre par mois. Aucun utilisateur sérieux d’une prévision ne s’en contente : un gestionnaire de barrage, un analyste de budget carbone et tous les centres météorologiques opérationnels — Météo-France comme le CEPMMT — demandent l’éventail des issues possibles et la probabilité de chacune. La machinerie tient en une petite modification de ce que nous avons déjà.
Un modèle apprend le quantile de l’issue lorsqu’il est entraîné sur la perte quantile (pinball loss)
qui facture τ par unité de sous-prévision et par unité de sur-prévision ; le minimiseur de son espérance est exactement le quantile τ. Prédisez et et vous obtenez un intervalle de prédiction à 90 % ; vous donne en prime une prévision ponctuelle médiane.
Nous dotons le LSTM de la section 3.4 d’une tête quantile : au lieu de 48 nombres, il en produit — un par pas d’horizon et par niveau de quantile — entraînée sur la perte quantile moyennée sur les niveaux. Mêmes fenêtres, même ancrage, mêmes époques. Cela nous achète une affirmation vérifiable, et c’est la vérification qui compte :
- Couverture empirique : sur les 48 mois de test, à quelle fréquence l’observation tombe-t-elle réellement dans l’intervalle à 90 % ? Une prévision calibrée couvre environ 90 %. Moins signifie que le modèle est trop confiant ; beaucoup plus signifie que les intervalles sont trop larges pour éclairer une décision. C’est la même discipline de calibration que les courbes de fiabilité des chapitres 3 et 4.5, appliquée aux intervalles.
- Un pli à défroisser : des sorties de quantiles entraînées indépendamment peuvent se croiser ( pour tel mois). Trier selon l’axe des quantiles au moment de la prédiction est la réparation standard en une ligne, et nous l’appliquons ci-dessous.
QUANTILES = np.round(np.arange(0.05, 0.951, 0.05), 2) # 0.05, 0.10, ..., 0.95
N_Q = len(QUANTILES)
q_torch = torch.tensor(QUANTILES, dtype=torch.float32, device=device).view(1, 1, -1)
class QuantileLSTMForecaster(nn.Module):
"""The 3.4 LSTM with a (horizon x n_quantiles) head."""
def __init__(self, horizon, n_q, hidden=32):
super().__init__()
self.horizon, self.n_q = horizon, n_q
self.lstm = nn.LSTM(input_size=1, hidden_size=hidden, batch_first=True)
self.head = nn.Linear(hidden, horizon * n_q)
def forward(self, x):
_, (h, _) = self.lstm(x)
return self.head(h[-1]).view(-1, self.horizon, self.n_q)
def pinball_loss(pred, target):
"""Mean pinball loss over all quantile levels.
pred: (batch, horizon, n_q); target: (batch, horizon)."""
diff = target.unsqueeze(-1) - pred
return torch.mean(torch.maximum(q_torch * diff, (q_torch - 1.0) * diff))
torch.manual_seed(0)
qlstm_co2 = QuantileLSTMForecaster(H_CO2, N_Q).to(device)
ds_q = torch.utils.data.TensorDataset(torch.from_numpy(Xc).unsqueeze(-1),
torch.from_numpy(Yc))
dl_q = torch.utils.data.DataLoader(ds_q, batch_size=32, shuffle=True)
opt_q = torch.optim.Adam(qlstm_co2.parameters(), lr=1e-3)
for _ in range(40):
qlstm_co2.train()
for xb, yb in dl_q:
xb, yb = xb.to(device), yb.to(device)
opt_q.zero_grad()
pinball_loss(qlstm_co2(xb), yb).backward()
opt_q.step()
# Predict the test window, un-anchor, and sort along the quantile axis
# (independently trained quantile outputs can cross; sorting is the repair).
ctx = train_co2.to_numpy()[-CONTEXT_CO2:]
x_q = torch.from_numpy(((ctx - ctx[-1]) / scale_co2)
.astype(np.float32)).reshape(1, -1, 1).to(device)
qlstm_co2.eval()
with torch.no_grad():
q_pred = qlstm_co2(x_q).cpu().numpy()[0] # (48, 19), scaled units
q_pred = np.sort(q_pred, axis=1) * scale_co2 + ctx[-1] # back to ppm
lo, med, hi = q_pred[:, 0], q_pred[:, N_Q // 2], q_pred[:, -1] # q05, q50, q95
inside = (test_co2.values >= lo) & (test_co2.values <= hi)
print(f"empirical coverage of the 90% interval: {inside.mean():.1%} "
f"({inside.sum()}/{len(inside)} test months inside)")
print(f"interval width: {hi[0] - lo[0]:.2f} ppm at 1 month -> "
f"{hi[-1] - lo[-1]:.2f} ppm at 48 months")
fig, ax = plt.subplots(figsize=(9, 4))
recent = y.loc["2017":]
ax.plot(recent.index, recent.values, color=COLORS["observations"], lw=1.5,
label="observations")
ax.fill_between(test_co2.index, lo, hi, color=COLORS["LSTM"], alpha=0.25,
label="90% interval (q05-q95)")
ax.plot(test_co2.index, med, color=COLORS["LSTM"], lw=1.2, ls="--",
label="median (q50)")
ax.axvline(test_co2.index[0], color="0.6", lw=0.8)
ax.set_xlabel("year")
ax.set_ylabel("CO$_2$ (ppm)")
ax.set_title("Quantile LSTM: 90% prediction interval on the test window")
ax.legend(fontsize=8)
fig.tight_layout()
plt.show()empirical coverage of the 90% interval: 91.7% (44/48 test months inside)
interval width: 2.33 ppm at 1 month -> 6.27 ppm at 48 months
Disons le chiffre sans détour : dans cette exécution, l’intervalle à 90 % attrape 44 des 48 mois de test — 91,7 % de couverture empirique pour une couverture nominale de 90 %. C’est un intervalle calibré, et rien ne le garantissait : la régression quantile vise la calibration mais rien ne l’impose hors échantillon, et c’est pourquoi on calcule la couverture au lieu de l’affirmer. Deux réserves honnêtes. D’abord, 48 mois consécutifs issus d’une origine de prévision unique sont fortement corrélés : cette estimation de couverture est donc elle-même bruitée — une autre fenêtre de test pourrait tout aussi bien renvoyer 80 % ou 100 %. Ensuite, regardez les largeurs affichées : l’intervalle passe d’environ 2,3 ppm à un mois d’échéance à environ 6,3 ppm à quatre ans. Le modèle a appris des fenêtres d’entraînement que l’incertitude s’accumule avec l’échéance ; personne ne le lui a dit.
Le CRPS. La couverture vérifie un intervalle à un niveau. Le CRPS (continuous ranked probability score, score de probabilité classé continu) — la règle de notation standard de la vérification des prévisions, le nombre que rapportent les systèmes d’ensemble atmosphériques — note la distribution prédictive tout entière face à l’issue :
Il a les unités des données (ici des ppm), et pour une prévision réduite à un point unique il se ramène exactement à l’erreur absolue — le CRPS moyen est donc directement comparable à la MAE d’une prévision ponctuelle, et il récompense d’un même geste la calibration et la finesse. La seconde identité est aussi l’estimateur : moyenner la perte quantile sur une grille dense de quantiles et doubler. Nous prédisons déjà 19 quantiles : l’implémentation tient en quelques lignes.
def crps_from_quantiles(y_true, q_pred, quantiles):
"""Quantile-based CRPS estimator: 2 x the pinball loss averaged over a
dense quantile grid. Exact in the limit of a dense grid; with 19 levels
it slightly undervalues the tails beyond q05/q95."""
yv = np.asarray(y_true, dtype=float).reshape(-1, 1)
diff = yv - np.asarray(q_pred, dtype=float)
tau = np.asarray(quantiles, dtype=float).reshape(1, -1)
return 2.0 * np.maximum(tau * diff, (tau - 1.0) * diff).mean(axis=1)
crps_q = crps_from_quantiles(test_co2.values, q_pred, QUANTILES)
print(f"mean CRPS, quantile LSTM : {crps_q.mean():.3f} ppm")
print(f"MAE of its median (q50) : {np.mean(np.abs(test_co2.values - med)):.3f} ppm")
print("point forecasts (CRPS = MAE):")
for name in ["LSTM", "SARIMA"]:
print(f" {name:12s}: {results_co2[name]['MAE']:.3f} ppm")mean CRPS, quantile LSTM : 0.766 ppm
MAE of its median (q50) : 1.051 ppm
point forecasts (CRPS = MAE):
LSTM : 0.886 ppm
SARIMA : 0.281 ppm
Lisez les trois nombres ensemble. Le CRPS moyen du LSTM quantile (0,766 ppm) bat la MAE de sa propre prévision médiane (1,051 ppm) et celle du LSTM à prévision ponctuelle de la section 3.4 (0,886 ppm) : étaler la probabilité sur la plage où l’issue peut tomber rapporte un vrai score, même quand le centre est imparfait. Cette comparaison est, en une ligne de sortie, tout l’argument de la prévision probabiliste. Mais la prévision ponctuelle de SARIMA l’emporte encore avec 0,281 ppm — un point bien placé bat une distribution honnêtement large. La calibration n’excuse pas un centre biaisé ou diffus ; le mot d’ordre de la communauté de la vérification est finesse sous contrainte de calibration, et le concurrent naturel suivant serait la propre distribution de prévision de SARIMA (get_forecast(...).conf_int() en donne les intervalles — essayez de les noter).
3.8 Erreur en fonction de l’échéance : l’horizon de skill¶
Une unique MAE à 48 mois écrase toute une dimension de la qualité d’une prévision : la vitesse à laquelle le skill (score de compétence) décroît avec l’échéance. La vérification atmosphérique ne rapporte jamais un seul nombre ; elle rapporte l’erreur en fonction de l’échéance, parce qu’un modèle excellent à 1 mois et inutile à 12 sert d’autres décisions qu’un modèle médiocre partout.
Nous la mesurons par une évaluation à origine glissante, la façon dont les systèmes opérationnels sont notés. Les modèles restent figés exactement tels qu’entraînés sur les données jusqu’à 2019-08 — aucun réentraînement. L’origine de prévision glisse ensuite mois par mois dans la fenêtre de test ; à chaque origine, le modèle voit l’historique observé jusqu’à ce mois (un prévisionniste opérationnel dispose toujours des dernières observations) et émet des prévisions aux échéances 1 à 12. En moyennant l’erreur absolue à chaque échéance sur toutes les origines, on obtient une courbe d’erreur en fonction de l’échéance — avec 37 origines par échéance au lieu de l’unique échantillon par échéance que fournit une seule prévision à 48 mois.
Le vocabulaire à retenir : l’horizon de skill (skill horizon) est l’échéance au-delà de laquelle un modèle cesse de battre la prévision de référence (ici la persistance). Au-delà, le modèle n’apporte rien de plus que « pas de changement ». Les modèles météorologiques atteignent leur horizon de skill vers deux semaines ; la question est de savoir où ces modèles-ci atteignent le leur sur le CO2.
LEADS = 12
LAGS_CO2, ROLL_CO2 = [1, 2, 3, 6, 12, 24], 12
y_vals = y.to_numpy()
first_test = len(train_co2)
origins = np.arange(first_test - 1, len(y) - LEADS) # index of the last observed month
# One LightGBM fit on the training series only, reused frozen at every origin.
X_tr_lead, y_tr_lead = build_training_table(train_co2, LAGS_CO2, ROLL_CO2)
lgbm_frozen = LGBMRegressor(n_estimators=300, learning_rate=0.05, num_leaves=15,
random_state=0, verbose=-1).fit(X_tr_lead, y_tr_lead)
abs_err = {name: np.zeros((len(origins), LEADS))
for name in ["naive", "seasonal naive", "LightGBM", "LSTM"]}
for r, j in enumerate(origins):
hist = y_vals[:j + 1] # everything observed at the origin
target = y_vals[j + 1:j + 1 + LEADS]
abs_err["naive"][r] = np.abs(target - hist[-1])
abs_err["seasonal naive"][r] = np.abs(target - y_vals[j + 1 - 12:j + 1 + LEADS - 12])
h = list(hist) # LightGBM: recursive from observed history
for step in range(LEADS):
row = pd.DataFrame([lag_feature_row(np.asarray(h), y.index[j + 1 + step],
LAGS_CO2, ROLL_CO2)])
h.append(h[-1] + float(lgbm_frozen.predict(row)[0]))
abs_err["LightGBM"][r] = np.abs(target - np.asarray(h[-LEADS:]))
abs_err["LSTM"][r] = np.abs(target - nn_forecast(lstm_co2, hist[-CONTEXT_CO2:],
scale_co2)[:LEADS])
leads = np.arange(1, LEADS + 1)
lead_mae = pd.DataFrame({name: e.mean(axis=0) for name, e in abs_err.items()},
index=leads)
naive_train_mae_co2 = np.mean(np.abs(np.diff(train_co2.values)))
fig, ax = plt.subplots(figsize=(7, 3.5))
for name in lead_mae.columns:
ax.plot(leads, lead_mae[name], color=COLORS[name], lw=1.5, marker="o",
ms=3, label=name)
sec = ax.secondary_yaxis("right", functions=(lambda v: v / naive_train_mae_co2,
lambda v: v * naive_train_mae_co2))
sec.set_ylabel("MASE")
ax.set_xlabel("lead time (months)")
ax.set_ylabel("MAE (ppm)")
ax.set_title(f"Error vs lead time, {len(origins)} rolling origins in the test window")
ax.legend(fontsize=8)
fig.tight_layout()
plt.show()
for name in ["seasonal naive", "LightGBM", "LSTM"]:
wins = leads[lead_mae[name].values < lead_mae["naive"].values]
print(f"{name:15s} beats persistence at leads: {wins.tolist()}")
seasonal naive beats persistence at leads: [3, 4, 5, 6, 7, 8, 9]
LightGBM beats persistence at leads: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
LSTM beats persistence at leads: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
La courbe et le contrôle affiché répondent à la question de l’horizon de skill, et sur le CO2 la réponse est l’inverse de celle de la météorologie. La persistance se dégrade régulièrement avec l’échéance — la tendance et le cycle saisonnier éloignent la série toujours plus de « pas de changement » — tandis que l’erreur des modèles entraînés croît bien plus lentement. LightGBM et le LSTM battent la persistance à toutes les échéances de 1 à 12 mois : leur horizon de skill se situe au-delà d’un an, plus loin que la plus longue échéance mesurée. Les modèles météorologiques perdent contre la climatologie en deux semaines environ parce que le chaos atmosphérique détruit l’information de l’état initial ; le CO2 mensuel conserve son skill parce que la structure qui porte la prévision (tendance plus cycle annuel) ne cesse jamais de s’appliquer. L’horizon de skill est une propriété de la série autant que du modèle.
Les prévisions de référence restent instructives. Le naïf saisonnier ne bat la persistance qu’aux échéances intermédiaires (3 à 9 dans cette exécution) : aux échéances 1–2, la valeur du mois précédent est simplement très proche de la vérité, et près de l’échéance 12 la persistance retombe dans la même saison, comblant l’essentiel de l’écart. Le choix même de la référence dépend donc de l’échéance — une raison de plus pour qu’un résumé en un seul nombre du skill d’une prévision cache plus qu’il ne montre. Pour une série où l’horizon de skill arrive presque tout de suite, appliquez cette même évaluation à origine glissante à la série de Jakobshavn de la section 4.
3.9 Le classement est-il stable ? La dispersion sur les graines¶
La section 4 avertira que les classements issus de séries courtes sont instables d’une graine aléatoire à l’autre. Le même aveu vaut dès ici, alors quantifions-le au lieu de le laisser à l’état de réserve. Tous les nombres du tableau 3.6 proviennent d’une seule exécution d’entraînement ; une entrée de tableau est un tirage dans une distribution sur les graines, et rapporter un tirage comme s’il s’agissait de la moyenne est le même péché que rapporter une prévision ponctuelle sans intervalle.
La sensibilité à la graine diffère selon la classe de modèle, et cela fait partie de la leçon :
- naïf, naïf saisonnier, SARIMA sont déterministes : pas de graine, pas de dispersion.
- LightGBM accepte un
random_state, mais tel qu’il est configuré ici (sans sous-échantillonnage de lignes ni de caractéristiques), l’algorithme est lui aussi déterministe — la graine est un bouton inerte tant que vous n’activez pas le bagging. Nous exécutons 5 graines pour le démontrer, pas pour estimer une dispersion. - Le LSTM et le transformer sont stochastiques deux fois : par l’initialisation des poids et par le brassage des lots (plus le dropout, pour le transformer). Ce sont ces entrées dont il faut se méfier de la position dans le tableau à partir d’une seule exécution. Ce sont aussi les entrées lentes : nous employons donc 3 graines chacune plutôt que 5, et nous le disons.
def mase_of(pred):
return forecast_metrics(test_co2.values, pred, train_co2.values)["MASE"]
seed_mase = {}
seed_mase["LightGBM"] = [
mase_of(lgbm_recursive_forecast(train_co2, test_co2.index,
lags=[1, 2, 3, 6, 12, 24], roll=12,
random_state=s)[0])
for s in range(5)
]
def deep_seed_mase(model_cls, seed):
torch.manual_seed(seed)
m = model_cls(horizon=H_CO2)
train_forecaster(m, Xc, Yc, epochs=40)
return mase_of(nn_forecast(m, train_co2.to_numpy()[-CONTEXT_CO2:], scale_co2))
# Seed 0 is the run already in the 3.6 table; two more per deep model. Three
# seeds, not five: these are the slow entries, and the spread they show is
# already unmistakable.
seed_mase["LSTM"] = ([results_co2["LSTM"]["MASE"]]
+ [deep_seed_mase(LSTMForecaster, s) for s in (1, 2)])
seed_mase["Transformer"] = ([results_co2["Transformer"]["MASE"]]
+ [deep_seed_mase(TransformerForecaster, s) for s in (1, 2)])
spread_table = pd.DataFrame({
name: {"MASE mean": np.mean(v), "MASE std": np.std(v),
"MASE min": np.min(v), "MASE max": np.max(v), "seeds": len(v)}
for name, v in seed_mase.items()}).T.round(3)
spread_tableLe tableau de dispersion affine le classement de la section 3.6 plutôt qu’il ne le renverse — et il montre exactement où loge le hasard :
- LightGBM affiche une MASE identique pour les cinq graines (écart-type 0,000), comme annoncé : sans sous-échantillonnage de lignes ni de caractéristiques, l’algorithme est déterministe, et « nous avons varié la graine » n’est une preuve de robustesse que lorsque l’algorithme consomme effectivement la graine.
- Le LSTM couvre 0,808–0,856 (moyenne 0,827, écart-type 0,021) ; le transformer couvre 0,327–0,403 (moyenne 0,356, écart-type 0,034) — un dixième de sa propre moyenne, du seul fait de l’initialisation, de l’ordre des lots et du dropout.
- Les étendues ne se recouvrent ni entre elles, ni avec le 0,256 déterministe de SARIMA, ni avec le 0,939 de LightGBM : sur cette exécution, l’ordre du concours survit donc aux graines. C’est une conclusion que nous pouvons désormais énoncer avec des preuves plutôt qu’avec de l’espoir — et ce n’est pas une loi générale. Sur la série de glace ci-dessous, les écarts entre modèles sont de la même taille que des dispersions comme celles-ci, raison pour laquelle son classement ne doit pas être cru sur la foi d’une exécution unique. Rapportez moyenne et dispersion ; laissez le lecteur voir si les écarts battent le bruit.
4. La série de glace : une épreuve plus dure¶
Nous ne faisons passer la série de glace que par les références, ARIMA, LightGBM et le LSTM. Moins de modèles, et c’est délibéré : la série est courte (environ 350 échantillons d’entraînement), échantillonnée irrégulièrement et non dominée par un cycle saisonnier net, si bien qu’un SARIMA saisonnier n’a pas grand-chose à saisir et qu’un transformer dispose de trop peu de données pour justifier ses paramètres. Les ajouter gonflerait le tableau, pas la compréhension.
4.1 Références et ARIMA¶
L’échantillonnage donne en moyenne environ 46 échantillons par an : le « naïf saisonnier » duplique donc les 46 derniers échantillons d’entraînement comme copie approximative de l’an dernier ; avec un échantillonnage irrégulier, cet alignement n’est qu’approximatif, raison de plus pour que la machinerie saisonnière aide moins ici. Pour l’entrée classique, nous employons un simple ARIMA(1,1,1), sans termes saisonniers.
from statsmodels.tsa.arima.model import ARIMA
H_ICE = len(test_ice)
m_ice = 46 # approx. one year (~46 samples/yr on average)
naive_ice = np.full(H_ICE, train_ice.iloc[-1])
seasonal_ice = np.tile(train_ice.values[-m_ice:],
int(np.ceil(H_ICE / m_ice)))[:H_ICE]
register(results_ice, preds_ice, "naive", naive_ice, test_ice, train_ice)
register(results_ice, preds_ice, "seasonal naive", seasonal_ice, test_ice, train_ice)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
arima_ice = ARIMA(train_ice.values, order=(1, 1, 1)).fit()
register(results_ice, preds_ice, "ARIMA", arima_ice.forecast(H_ICE), test_ice, train_ice)naive -> {'MSE': np.float64(34556.717), 'RMSE': np.float64(185.894), 'MAE': np.float64(168.382), 'MAPE (%)': np.float64(9.589), 'MASE': np.float64(6.249)}
seasonal naive -> {'MSE': np.float64(23662.811), 'RMSE': np.float64(153.827), 'MAE': np.float64(129.081), 'MAPE (%)': np.float64(7.237), 'MASE': np.float64(4.79)}
ARIMA -> {'MSE': np.float64(32800.223), 'RMSE': np.float64(181.108), 'MAE': np.float64(163.157), 'MAPE (%)': np.float64(9.297), 'MASE': np.float64(6.055)}
4.2 LightGBM et le LSTM, réutilisés¶
Mêmes fonctions utilitaires, nouveaux jeux de retards : des retards jusqu’à 46 échantillons (environ un an) et une moyenne glissante courte. Le LSTM emploie un contexte d’environ un an (46 échantillons) et prévoit directement toute la fenêtre de test, comme précédemment.
lgbm_pred_ice, _ = lgbm_recursive_forecast(train_ice, test_ice.index,
lags=[1, 2, 3, 6, 23, 46], roll=6)
register(results_ice, preds_ice, "LightGBM", lgbm_pred_ice, test_ice, train_ice)LightGBM -> {'MSE': np.float64(13027.719), 'RMSE': np.float64(114.139), 'MAE': np.float64(94.373), 'MAPE (%)': np.float64(5.338), 'MASE': np.float64(3.502)}
CONTEXT_ICE = 46
Xi, Yi = make_windows(train_ice.to_numpy(), CONTEXT_ICE, H_ICE)
scale_ice = float(Xi.std())
Xi, Yi = Xi / scale_ice, Yi / scale_ice
print(f"{len(Xi)} training windows, scale = {scale_ice:.0f} m/yr")
torch.manual_seed(0)
lstm_ice = LSTMForecaster(horizon=H_ICE)
_ = train_forecaster(lstm_ice, Xi, Yi, epochs=40)
lstm_pred_ice = nn_forecast(lstm_ice, train_ice.to_numpy()[-CONTEXT_ICE:], scale_ice)
register(results_ice, preds_ice, "LSTM", lstm_pred_ice, test_ice, train_ice)213 training windows, scale = 83 m/yr
LSTM -> {'MSE': np.float64(52166.004), 'RMSE': np.float64(228.399), 'MAE': np.float64(191.047), 'MAPE (%)': np.float64(10.931), 'MASE': np.float64(7.09)}
ice_table = pd.DataFrame(results_ice).T.sort_values("MASE").round(3)
ice_tablefig, ax = plt.subplots(figsize=(9, 4))
recent_ice = ice.loc["2020":]
ax.plot(recent_ice.index, recent_ice.values, color=COLORS["observations"],
lw=1.5, label="observations")
for name, pred in preds_ice.items():
ax.plot(test_ice.index, pred, color=COLORS[name], lw=1.2, ls="--", label=name)
ax.axvline(test_ice.index[0], color="0.6", lw=0.8)
ax.set_xlabel("year")
ax.set_ylabel("surface speed (m/yr)")
ax.set_title("Ice test window: all forecasts vs observations")
ax.legend(fontsize=8, ncol=2)
fig.tight_layout()
plt.show()
Notez à quel point les valeurs de MASE sont plus grandes que pour le CO2, et ce pour tous les modèles. Le jeu naïf à un pas (le dénominateur de la MASE) est facile sur une série de vitesse à variation lente, mais la fenêtre de test pluriannuelle contient des changements de comportement du glacier qu’aucun motif de la fenêtre d’entraînement n’annonce. Quand l’avenir n’est pas écrit dans le passé, plus de modèle ne veut pas dire plus de prévision. Attendez-vous aussi à ce que le classement entre modèles entraînés soit instable ici : avec ~340 échantillons d’entraînement, de petits choix (longueur du contexte, jeu de retards, graine aléatoire) font monter et descendre les modèles dans le tableau — nous avons mesuré exactement cet effet sur le concours CO2 à la section 3.9, et c’est pire sur une série aussi courte. Rapportez-le honnêtement plutôt que de partir en quête des réglages qui font le meilleur effet.
5. Les événements rares cassent les métriques globales¶
Le problème de données emblématique des géosciences, c’est un champ qui varie doucement pendant des années puis fait quelque chose d’extrême pendant trois jours : les crues sur un enregistrement de débit, les éruptions sur un enregistrement de trémor, les sauts sur un enregistrement de déformation. Toutes les métriques de ce carnet moyennent jusqu’ici sur tous les échantillons de test, et les échantillons calmes sont environ douze fois plus nombreux que les échantillons d’événement — un modèle peut donc afficher une excellente MAE tout en ratant chaque événement qui compte.
Pour le montrer proprement, il nous faut une vérité terrain sur les événements : nous construisons donc la série nous-mêmes avec mlgeo_synth.inject_rare_events, soit huit ans d’un fond saisonnier journalier, ponctué d’environ 6 événements par an durant 2 à 6 jours, dont les amplitudes sont tirées d’une loi de Pareto généralisée (tail="gpd", paramètre de forme ). C’est une queue en loi de puissance — le régime des valeurs extrêmes des enregistrements réels de crues et de surcotes, où la variance est infinie et où le plus grand événement d’une fenêtre future à venir est couramment plusieurs fois plus grand que le plus grand observé jusqu’ici. Le générateur renvoie une colonne amplitude : la vraie taille tirée pour chaque événement, exactement la vérité terrain qu’une métrique globale ne voit jamais.
Le prévisionniste reçoit délibérément tous les avantages : un LightGBM à un pas d’avance (la valeur de la veille est toujours dans le jeu de caractéristiques), entraîné sur six ans, testé sur les deux dernières années. Puis, au lieu d’une MAE unique, nous stratifions l’erreur de test par l’amplitude de vérité terrain.
n_days_ev = int(8 * 365.25)
days_ev = pd.date_range("2015-01-01", periods=n_days_ev, freq="D")
rng_ev = np.random.default_rng(0)
background_ev = (10.0 + 3.0 * np.sin(2 * np.pi * np.arange(n_days_ev) / 365.25)
+ 0.3 * rng_ev.standard_normal(n_days_ev))
events = mlgeo_synth.inject_rare_events(background_ev, rate_per_year=6.0,
duration_days=(2, 6), shape="spike",
tail="gpd", seed=7)
series_ev = pd.Series(events["value"].to_numpy(), index=days_ev)
# One-step-ahead LightGBM on lag features — the easiest possible forecasting
# task: the model always has yesterday's observation in hand.
LAGS_EV, ROLL_EV = [1, 2, 3, 5, 7, 14, 365], 7
X_ev, dy_ev = build_training_table(series_ev, LAGS_EV, ROLL_EV)
start_ev = max(max(LAGS_EV), ROLL_EV)
dates_ev = series_ev.index[start_ev:]
cut = dates_ev.searchsorted(pd.Timestamp("2021-01-01")) # train < 2021 <= test
lgbm_ev = LGBMRegressor(n_estimators=300, learning_rate=0.05, num_leaves=15,
random_state=0, verbose=-1).fit(X_ev.iloc[:cut], dy_ev[:cut])
pred_ev = (series_ev.to_numpy()[start_ev - 1:-1][cut:]
+ lgbm_ev.predict(X_ev.iloc[cut:])) # previous value + predicted change
true_ev = series_ev.to_numpy()[start_ev:][cut:]
err_ev = np.abs(true_ev - pred_ev)
# Stratify by the generator's ground-truth amplitude column.
evt = events.iloc[start_ev:].iloc[cut:] # aligned truth columns
quiet = evt["event"].to_numpy() == 0
amp = evt["amplitude"].to_numpy()
strata = {
"all test samples": np.ones_like(quiet),
"quiet (no event)": quiet,
"event, amplitude < 10": ~quiet & (amp < 10),
"event, amplitude 10-25": ~quiet & (amp >= 10) & (amp < 25),
"event, amplitude >= 25": ~quiet & (amp >= 25),
}
tail_table = pd.DataFrame({name: {"n samples": int(m.sum()),
"MAE": err_ev[m].mean()}
for name, m in strata.items()}).T
tail_table["n samples"] = tail_table["n samples"].astype(int)
print(tail_table.round(3))
print(f"\nlargest single miss: {err_ev.max():.1f} "
f"(event amplitude {amp[np.argmax(err_ev)]:.1f})")
# Left: the test window. Right: per-event peak error vs true amplitude.
per_event = (pd.DataFrame({"event_id": evt["event_id"].to_numpy(),
"amp": amp, "err": err_ev})
.query("event_id >= 0").groupby("event_id")
.agg(amp=("amp", "max"), peak_err=("err", "max")))
fig, axes = plt.subplots(1, 2, figsize=(10, 3.5),
gridspec_kw={"width_ratios": [2, 1]})
axes[0].plot(dates_ev[cut:], true_ev, color=COLORS["observations"], lw=0.8,
label="observed")
axes[0].plot(dates_ev[cut:], pred_ev, color=COLORS["LightGBM"], lw=0.8,
alpha=0.8, label="one-step forecast")
axes[0].set_ylabel("value")
axes[0].set_title("Test window (events ride on the seasonal background)")
axes[0].legend(fontsize=8)
axes[1].loglog(per_event["amp"], per_event["peak_err"], "o", ms=5,
color=COLORS["LightGBM"])
lim = [4, per_event["amp"].max() * 1.5]
axes[1].plot(lim, lim, color="0.6", lw=0.8, ls="--", label="error = amplitude")
axes[1].set_xlabel("true event amplitude")
axes[1].set_ylabel("peak |error| in event")
axes[1].set_title("Misses scale with amplitude")
axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show() n samples MAE
all test samples 730 0.969
quiet (no event) 674 0.621
event, amplitude < 10 43 2.458
event, amplitude 10-25 7 10.577
event, amplitude >= 25 6 18.274
largest single miss: 50.8 (event amplitude 50.2)
Une lecture naïve de la première ligne dit que le modèle est excellent : une MAE de 0,97 sur une série dont la seule oscillation saisonnière fait 6 unités. Les lignes stratifiées racontent la vraie histoire. Les jours calmes obtiennent 0,62 — le plancher de bruit, le meilleur que puisse faire une prévision à un pas. Les petits événements (amplitude inférieure à 10) coûtent déjà 2,5. Les sept échantillons des événements moyens coûtent environ 10,6, les six échantillons des plus grands événements 18,3, et le plus grand événement de tous — amplitude 50,2 — est raté de 50,8 : le modèle a prédit le fond saisonnier et la crue a eu lieu quand même. Le panneau de droite rend la chose systématique : l’erreur maximale par événement suit la vraie amplitude le long de la droite erreur = amplitude. Ce que le modèle n’a jamais vu, il le rate en entier.
Les deux causes sont structurelles, ce ne sont pas des bogues. Les déclenchements d’événements ne sont pas prévisibles à partir de caractéristiques de retard — aucune valeur de la veille ne vous dit qu’un pic commence demain — et la queue de Pareto généralisée () garantit que la fenêtre de test contient des amplitudes au-delà de tout exemple d’entraînement : le modèle extrapole donc précisément sur les échantillons qui comptent le plus. Pendant ce temps, la métrique globale masque tout cela : les échantillons d’événement représentent 56 des 730 jours de test (environ un sur treize), si bien que la MAE globale (0,97) reste à un facteur 1,6 du plancher des jours calmes, alors que la strate supérieure est trente fois pire. Un classement noté sur la MAE globale féliciterait ce modèle.
La réparation est dans la métrique, pas dans le modèle : si votre prévision existe pour les événements — crues, éruptions, sauts —, notez-les séparément. Des tableaux stratifiés comme celui-ci, un rappel par fenêtre d’événement ou des pertes pondérées par l’amplitude font tous l’affaire ; moyenner sur les jours calmes, non. C’est la même discipline d’évaluation équitable que la piste cachée du classement ci-dessous : ce que vous mesurez est ce que vous optimisez, alors mesurez ce qui compte.
6. Exercices¶
Exercice 1. Dans l’entrée LightGBM du CO2, retirez lag_12 et lag_24 de la liste des retards et relancez. Qu’advient-il de la MASE, et pourquoi ?
Solution
La MASE se dégrade. Les retards de 12 et 24 mois sont la façon dont le modèle voit le cycle annuel ; sans eux, il doit reconstruire la saisonnalité à partir de la seule caractéristique de mois calendaire et de retards courts, et la prévision récursive dérive lentement hors de la forme saisonnière. Exécutez
lgbm_recursive_forecast(train_co2, test_co2.index, lags=[1, 2, 3, 6], roll=12)
et comparez les métriques avec forecast_metrics.
Exercice 2. Les modèles profonds prévoient les 48 mois d’un coup (multi-pas direct) ; LightGBM prédit un mois à la fois et réinjecte ses prédictions (récursif). Citez un avantage de chaque stratégie.
Solution
Direct : pas de rétroaction d’erreur, donc un mauvais pas précoce ne peut pas contaminer les pas suivants, et le modèle est entraîné exactement sur l’horizon sur lequel il est noté. Récursif : un seul modèle sert n’importe quel horizon, et l’entraînement à un pas utilise bien plus d’exemples par paramètre ; mais les erreurs se composent le long de l’horizon. Des hybrides (un modèle direct par pas, ou les schémas dits « rectify ») existent précisément pour cet arbitrage.
Exercice 3. Relancez la section 5 deux fois : une fois avec tail="uniform" (le défaut du générateur : amplitudes uniformes sur 5–20) et une fois en gardant tail="gpd" mais en portant rate_per_year à 24. Lequel des deux changements fait le plus pour la strate des grands événements, et que vous apprend cela sur la différence entre rare et à queue lourde ?
Solution
La queue uniforme aide davantage. Avec des amplitudes plafonnées à 20, la fenêtre d’entraînement a déjà montré au modèle des événements aussi grands que tous ceux sur lesquels il sera noté — le problème est de l’interpolation, et davantage d’exemples (taux plus élevé) la peaufinent. Avec une queue de Pareto généralisée, le plus grand événement futur est hors de portée de tout exemple d’entraînement presque par construction : aucun taux d’événements réaliste ne répare la strate supérieure, car on demande au modèle d’extrapoler. La rareté limite le nombre d’exemples que vous obtenez ; une queue lourde garantit que l’examen contient une magnitude que le manuel n’a jamais montrée. (L’exercice hors distribution du chapitre 4.5 est la même leçon dans l’espace des caractéristiques.)
7. Le classement de la classe : deux pistes¶
Les 12 mois retirés à la section 1.2 reviennent maintenant en jeu — mais d’abord, une déclaration honnête sur ce qu’un classement public peut et ne peut pas mesurer.
Piste A — CO2 (diagnostic).
- Produisez une prévision de CO2 à 12 mois pour 2023-09 à 2024-08, les mois réservés, avec n’importe quel modèle entraîné sur les données allant jusqu’à 2023-08 seulement. Sur l’honneur — voyez l’encadré ci-dessus pour comprendre pourquoi tricher ici ne fait que ruiner le diagnostic.
- Enregistrez-la dans
results/forecast_<uwnetid>.csvavec exactement deux colonnes :date(au formatYYYY-MM) etco2_ppm. - Soumettez-la par une pull request sur le dépôt du livre du cours, https://
github .com /geo -smart /mlgeo -book. L’intégration continue note chaque PR et publie le classement.
La note est la MASE face aux mois réservés, mise à l’échelle par la MAE naïve à un pas de la série d’entraînement (la même convention que forecast_metrics à la section 2.2 — le dénominateur vient de l’enregistrement antérieur aux données réservées, pas des données réservées elles-mêmes).
Piste B — la série synthétique cachée (notée). Le fichier leaderboard/synth_forecast_history.csv du dépôt du cours contient une série journalière de déplacement produite par mlgeo_synth.gnss_series — tendance séculaire, cycles annuel et semi-annuel, et bruit coloré, à l’image des enregistrements GNSS du chapitre 2. Elle est générée à partir d’une graine privée et régénérée chaque année. Prévoyez les 90 jours suivant la dernière date de l’historique et enregistrez results/forecast_hidden_<uwnetid>.csv avec les colonnes date (YYYY-MM-DD) et disp_mm, soumis dans la même pull request. La série de vérité n’est jamais versionnée ; le classement public indique votre soumission comme reçue, et c’est l’exécution de l’enseignante qui produit le tableau noté. Même métrique : la MASE, mise à l’échelle par l’erreur naïve à un pas de l’historique.
Les cellules ci-dessous écrivent un fichier d’exemple pour chaque piste à l’aide du modèle LightGBM récursif. Remplacez example par votre UW NetID et les prévisions par les vôtres.
uwnetid = "example" # <-- replace with your UW NetID
# Track A: refit on ALL pre-holdout data (through 2023-08) and forecast the
# 12 holdout months.
final_pred, _ = lgbm_recursive_forecast(y, holdout_months,
lags=[1, 2, 3, 6, 12, 24], roll=12)
submission = pd.DataFrame({"date": holdout_months.strftime("%Y-%m"),
"co2_ppm": np.round(final_pred, 2)})
os.makedirs("results", exist_ok=True)
out_path = f"results/forecast_{uwnetid}.csv"
submission.to_csv(out_path, index=False)
print("wrote", out_path)
print(submission.to_string(index=False))wrote results/forecast_example.csv
date co2_ppm
2023-09 418.23
2023-10 418.42
2023-11 419.86
2023-12 421.29
2024-01 422.44
2024-02 423.25
2024-03 424.07
2024-04 425.49
2024-05 426.23
2024-06 425.65
2024-07 423.90
2024-08 421.89
# Track B: forecast the hidden synthetic series 90 days past its history.
hist_path = Path("../../leaderboard/synth_forecast_history.csv")
if not hist_path.exists(): # running outside the repo checkout
hist_path = Path(pooch.retrieve(
"https://raw.githubusercontent.com/geo-smart/mlgeo-book/main/"
"leaderboard/synth_forecast_history.csv",
known_hash=None, fname="synth_forecast_history.csv"))
hidden_hist = pd.read_csv(hist_path, parse_dates=["date"])
hidden_series = pd.Series(hidden_hist["disp_mm"].to_numpy(),
index=pd.DatetimeIndex(hidden_hist["date"]))
horizon_dates = pd.date_range(hidden_series.index[-1] + pd.Timedelta(days=1),
periods=90, freq="D")
print(f"history: {len(hidden_series)} days ending {hidden_series.index[-1].date()}; "
f"forecast {horizon_dates[0].date()} .. {horizon_dates[-1].date()}")
hidden_pred, _ = lgbm_recursive_forecast(hidden_series, horizon_dates,
lags=[1, 2, 3, 7, 14, 365], roll=7)
hidden_sub = pd.DataFrame({"date": horizon_dates.strftime("%Y-%m-%d"),
"disp_mm": np.round(hidden_pred, 3)})
hidden_out = f"results/forecast_hidden_{uwnetid}.csv"
hidden_sub.to_csv(hidden_out, index=False)
print("wrote", hidden_out)
hidden_sub.head()history: 2832 days ending 2022-10-02; forecast 2022-10-03 .. 2022-12-31
wrote results/forecast_hidden_example.csv
Résumé¶
- Les références d’abord : les prévisions naïve et naïve saisonnière fixent la barre, et la MASE mesure chaque modèle face au jeu naïf à un pas.
- Un seul découpage temporel, une seule fonction de métriques, un seul tableau programmatique. Les découpages aléatoires et les nombres saisis à la main sont ce qui a fait déraper l’édition 2024 ; ne répétez ni l’un ni l’autre.
- Sur des séries longues et structurées (CO2), un modèle classique bien choisi (SARIMA) est difficile à battre, et les modèles profonds demandent du soin (ancrage, mise à l’échelle) rien que pour rivaliser. Sur des séries courtes et irrégulières (Jakobshavn), tous les modèles peinent, et le dire est le résultat correct.
- Le choix entre prévision multi-pas récursive et directe est un vrai choix de conception ; nous avons employé le récursif pour LightGBM et le direct pour le LSTM et le transformer.
- Une prévision ponctuelle est un début, pas une fin. Une tête quantile entraînée sur la perte quantile a produit un intervalle à 90 % dont nous avons mesuré la couverture empirique au lieu de la supposer, et le CRPS — la généralisation de la MAE propre à la vérification — a noté la distribution entière.
- L’erreur est une fonction de l’échéance ; tracez-la. L’horizon de skill — l’échéance à laquelle un modèle cesse de battre la persistance — est le résumé opérationnel d’un système de prévision.
- Un classement issu d’une seule graine est un seul tirage. Les entrées déterministes (références, SARIMA, cette configuration de LightGBM) ont une dispersion nulle ; les entrées profondes bougent assez d’une graine à l’autre pour échanger leurs places. Rapportez moyenne et dispersion, pas un tirage unique.
- Une MAE globale peut paraître excellente alors que tous les grands événements rares sont ratés : avec des amplitudes à queue lourde (Pareto généralisée), la strate supérieure d’événements portait des erreurs trente fois supérieures à la MAE des jours calmes, et le plus grand événement a été raté presque en entier. Stratifiez votre notation, ou votre métrique cachera précisément les échantillons qui comptent.
- Le classement comporte deux pistes parce que ses données réservées publiques sont des données publiques : le score CO2 est un diagnostic (et trivialement contournable — nous l’écrivons noir sur blanc), tandis que la série synthétique cachée, régénérée chaque année à partir d’une graine privée, porte le poids de la note.
Ensuite : soumettez vos deux prévisions pour le classement, puis passez au sujet du projet final (4.20). Pour les architectures de séquence elles-mêmes, revenez au carnet 4.4 ; pour la machinerie du boosting, au chapitre 3.