4.4 Modelos de secuencias: RNN, LSTM y atención - Machine Learning en las Geociencias

Las redes densas tratan cada característica de entrada como independiente. Las series de tiempo no lo son: el orden de las muestras lleva la señal. Este cuaderno construye cuatro modelos que leen una secuencia, en orden de sofisticación creciente: una red recurrente simple, una LSTM, una capa de autoatención de una sola cabeza escrita desde cero y un codificador transformer pequeño. Los cuatro resuelven la misma tarea de pronóstico sobre los mismos datos, así que podemos compararlos directamente.

Una neurona recurrente recibe una entrada y la salida que produjo en el paso de tiempo anterior. Como cada paso reutiliza la salida del paso previo, la red tiene memoria. Las celdas recurrentes simples tienen memoria corta, del orden de decenas de pasos; buena parte de este cuaderno trata de hacerlo mejor.

De Dive into Deep Learning: una RNN con estado oculto. En cada paso de tiempo la celda combina la entrada actual con el estado oculto del paso anterior.

🖥️ Diapositivas — Sesión 23 (lun 23 nov)

import time

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch

import mlgeo_synth

torch.manual_seed(42)
np.random.seed(42)

device = torch.device("cuda" if torch.cuda.is_available()
                      else "mps" if torch.backends.mps.is_available()
                      else "cpu")
print(f"Detected device: {device}")

# The models in this notebook are tiny, and recurrent layers launch one small
# operation per time step, so accelerator overhead dominates any speedup.
# CPU is usually fastest for models this size.
device = torch.device("cpu")
print(f"Using device: {device}")

Detected device: mps
Using device: cpu

1. Vocabulario de pronóstico¶

En el pronóstico de series de tiempo, la meta es predecir valores futuros a partir de datos históricos. Dos definiciones aparecen constantemente:

Ventana de contexto (también llamada ventana de historia o ventana de entrada): el tramo de datos pasados que el modelo lee antes de hacer una predicción. Su longitud fija cuánta historia puede usar el modelo para captar tendencias y estacionalidad.
Horizonte de pronóstico (también llamado longitud de predicción): el número de pasos de tiempo futuros que el modelo predice. Un horizonte de 30 significa que el modelo produce los siguientes 30 valores de una sola vez.

Usaremos una ventana de contexto de 90 días y un horizonte de pronóstico de 30 días.

2. Datos: una serie sintética de desplazamiento GNSS¶

Las estaciones GNSS registran la posición del suelo con precisión milimétrica, día tras día, durante décadas. Una serie de tiempo de desplazamiento en una estación mezcla varias señales físicas: movimiento tectónico constante, carga estacional de agua y nieve, saltos cosísmicos súbitos cuando ocurre un sismo y relajación postsísmica lenta después, todo sobre ruido instrumental y ambiental.

El paquete mlgeo_synth genera una serie así con componentes conocidos, de modo que podemos ver exactamente qué se les pide pronosticar a los modelos. Generamos 10 años de desplazamiento diario en milímetros, con un sismo colocado cerca de la mitad del registro.

eq_day = 1800  # earthquake day, near the middle of the 10-year record
gnss = mlgeo_synth.gnss_series(n_years=10.0, eq_day=eq_day, seed=42)
print(f"{len(gnss)} daily samples, columns: {list(gnss.columns)}")
gnss.head()

3652 daily samples, columns: ['date', 'disp_mm', 'trend_mm', 'seasonal_mm', 'eq_mm']

fig, axes = plt.subplots(2, 1, figsize=(8, 5), sharex=True)

axes[0].plot(gnss["date"], gnss["disp_mm"], color="#333333", lw=0.6,
             label="disp_mm (observed)")
axes[0].set_ylabel("displacement (mm)")
axes[0].legend(loc="upper left")
axes[0].grid(alpha=0.3)

axes[1].plot(gnss["date"], gnss["trend_mm"], color="#4477AA", label="trend_mm")
axes[1].plot(gnss["date"], gnss["seasonal_mm"], color="#CCBB44", label="seasonal_mm")
axes[1].plot(gnss["date"], gnss["eq_mm"], color="#EE6677", label="eq_mm")
axes[1].set_ylabel("component (mm)")
axes[1].set_xlabel("date")
axes[1].legend(loc="upper left")
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

El salto del sismo y la recuperación postsísmica curvada se montan sobre una tendencia constante de unos 12 mm/año y un ciclo estacional de unos pocos milímetros. La tarea de pronóstico: dados los últimos 90 días de desplazamiento, predecir los 30 días siguientes.

3. Pares supervisados y una división temporal¶

Cuatro reglas de preprocesamiento para el pronóstico de secuencias:

Estandarice la serie (reste la media, divida por la desviación estándar) para que la red entrene con valores cercanos a cero. Calcule las estadísticas solo con la porción de entrenamiento y luego aplíquelas en todas partes.
Deslice una ventana sobre la serie para construir pares (entrada, objetivo): 90 días de contexto como entrada, los 30 días siguientes como objetivo.
Divida en el tiempo, no al azar. Los primeros 8 años se vuelven datos de entrenamiento, los últimos 2 años validación. Ninguna ventana cruza la frontera de la división, y ningún dato se baraja a través de ella. Barajar ventanas dentro del conjunto de entrenamiento durante el entrenamiento está bien.
Ancle cada ventana. Reste el último valor del contexto tanto al contexto como al objetivo, de modo que el modelo pronostique el cambio relativo a la observación más reciente. Esta serie sube en tendencia durante 10 años, así que las ventanas de validación se ubican en niveles absolutos que el conjunto de entrenamiento nunca contuvo; sin el anclaje, todos los modelos tendrían que extrapolar fuera de su rango de entrenamiento y los que saturan (tanh, sigmoide) fallarían gravemente. Pruebe más adelante quitar el anclaje y observe crecer los errores de validación.

Una división aleatoria filtraría información: una ventana de «validación» podría solaparse casi por completo con dos ventanas de entrenamiento a ambos lados, y el puntaje de validación no diría nada sobre pronosticar tiempo genuinamente no visto.

WINDOW, HORIZON = 90, 30

series = gnss["disp_mm"].to_numpy(dtype=np.float32)
n_days = len(series)
split = int(0.8 * n_days)  # first 8 years train, last 2 years validation

mu = float(series[:split].mean())
sigma = float(series[:split].std())
z = (series - mu) / sigma
print(f"training-portion mean {mu:.1f} mm, std {sigma:.1f} mm")


def make_windows(z, start, stop, window=WINDOW, horizon=HORIZON):
    """Build anchored (window, horizon) pairs from z[start:stop]."""
    X, Y = [], []
    for i in range(start, stop - window - horizon + 1):
        X.append(z[i : i + window])
        Y.append(z[i + window : i + window + horizon])
    X = np.stack(X)[..., np.newaxis]      # (n_pairs, window, 1)
    Y = np.stack(Y)                       # (n_pairs, horizon)
    anchor = X[:, -1:, 0].copy()          # last context value, (n_pairs, 1)
    X = X - anchor[:, :, np.newaxis]      # context relative to its last day
    Y = Y - anchor                        # target relative to the same day
    return (torch.from_numpy(X), torch.from_numpy(Y),
            torch.from_numpy(anchor))


X_train, Y_train, anc_train = make_windows(z, 0, split)
X_val, Y_val, anc_val = make_windows(z, split, n_days)
print("train:", tuple(X_train.shape), "->", tuple(Y_train.shape))
print("val:  ", tuple(X_val.shape), "->", tuple(Y_val.shape))

training-portion mean 64.4 mm, std 46.6 mm
train: (2802, 90, 1) -> (2802, 30)
val:   (612, 90, 1) -> (612, 30)

i = 300  # one validation pair, shown in original units
a = anc_val[i].item()
ctx_mm = (X_val[i, :, 0].numpy() + a) * sigma + mu
tgt_mm = (Y_val[i].numpy() + a) * sigma + mu

plt.figure(figsize=(6, 3))
plt.plot(np.arange(WINDOW), ctx_mm, color="#4477AA", label="context (90 days)")
plt.plot(np.arange(WINDOW, WINDOW + HORIZON), tgt_mm, color="#EE6677",
         label="target (30 days)")
plt.axvline(WINDOW, color="gray", lw=0.8, ls="--")
plt.xlabel("day within window")
plt.ylabel("displacement (mm)")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

4. Un solo bucle de entrenamiento, una sola métrica¶

Cada modelo de abajo mapea un contexto de (90, 1) a 30 valores de pronóstico. Para mantener justa la comparación, todos comparten las mismas ventanas, la misma pérdida (MSE sobre valores estandarizados), la misma configuración del optimizador y la misma métrica de validación: error absoluto medio sobre el horizonte de 30 días, convertido de vuelta a milímetros.

El entrenamiento corre durante 20 épocas, suficiente para modelos así de pequeños. Aumente el número de épocas en su propia máquina si quiere una convergencia más ajustada.

results = {}


def count_params(model):
    return sum(p.numel() for p in model.parameters())


def val_mae_mm(model):
    """Mean absolute error on the validation windows, in millimeters."""
    model.eval()
    with torch.no_grad():
        pred = model(X_val.to(device)).cpu()
    return (pred - Y_val).abs().mean().item() * sigma


def train_model(model, name, n_epochs=20, lr=1e-3, batch_size=64):
    model = model.to(device)
    opt = torch.optim.Adam(model.parameters(), lr=lr)
    loss_fn = torch.nn.MSELoss()
    n = X_train.shape[0]
    history = {"train": [], "val": []}
    t0 = time.time()
    for epoch in range(n_epochs):
        model.train()
        perm = torch.randperm(n)  # shuffle within the training set only
        running = 0.0
        for j in range(0, n, batch_size):
            idx = perm[j : j + batch_size]
            xb, yb = X_train[idx].to(device), Y_train[idx].to(device)
            opt.zero_grad()
            loss = loss_fn(model(xb), yb)
            loss.backward()
            opt.step()
            running += loss.item() * len(idx)
        history["train"].append(running / n)
        model.eval()
        with torch.no_grad():
            val_loss = loss_fn(model(X_val.to(device)), Y_val.to(device)).item()
        history["val"].append(val_loss)
        if (epoch + 1) % 5 == 0:
            print(f"[{name}] epoch {epoch + 1:2d}  "
                  f"train MSE {history['train'][-1]:.4f}  val MSE {val_loss:.4f}")
    elapsed = time.time() - t0
    mae = val_mae_mm(model)
    results[name] = {"model": model, "history": history,
                     "n_params": count_params(model),
                     "train_time_s": elapsed, "val_mae_mm": mae}
    print(f"[{name}] {count_params(model):,} parameters, "
          f"trained in {elapsed:.1f} s, val MAE {mae:.2f} mm")

5. RNN simple¶

torch.nn.RNN implementa la celda recurrente simple. En cada paso de tiempo actualiza el estado oculto como $h_t = \tanh(W_x x_t + W_h h_{t-1} + b)$ : la entrada nueva mezclada con el estado oculto anterior. Leemos el estado oculto después del último de los 90 pasos y lo mapeamos a 30 valores de pronóstico con una capa lineal.

class VanillaRNN(torch.nn.Module):
    def __init__(self, hidden_size=32, horizon=HORIZON):
        super().__init__()
        self.rnn = torch.nn.RNN(1, hidden_size, batch_first=True)
        self.head = torch.nn.Linear(hidden_size, horizon)

    def forward(self, x):                # x: (batch, window, 1)
        out, _ = self.rnn(x)             # out: (batch, window, hidden)
        return self.head(out[:, -1, :])  # last hidden state -> 30 values


train_model(VanillaRNN(), "RNN")

[RNN] epoch  5  train MSE 0.0044  val MSE 0.0027

[RNN] epoch 10  train MSE 0.0042  val MSE 0.0025

[RNN] epoch 15  train MSE 0.0042  val MSE 0.0026

[RNN] epoch 20  train MSE 0.0042  val MSE 0.0025
[RNN] 2,110 parameters, trained in 4.5 s, val MAE 1.83 mm

hist = results["RNN"]["history"]
plt.figure(figsize=(6, 3))
plt.plot(hist["train"], color="#4477AA", label="train")
plt.plot(hist["val"], color="#EE6677", label="validation")
plt.yscale("log")
plt.xlabel("epoch")
plt.ylabel("MSE (standardized units)")
plt.title("Vanilla RNN learning curves")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

6. Por qué las RNN simples olvidan: gradientes desvanecidos y explosivos¶

Entrenar una red recurrente usa la retropropagación a través del tiempo (BPTT, por sus siglas en inglés): se desenrolla la red sobre los 90 pasos de tiempo y luego se empuja el gradiente de la pérdida hacia atrás a través de cada paso. El gradiente que llega al paso 1 es un producto de aproximadamente 90 matrices jacobianas, una por paso. Multiplique 90 números apenas menores que uno y el producto es casi cero; multiplique 90 números apenas mayores que uno y explota. Los productos de matrices se comportan igual: la multiplicación repetida se encoge o se dispara exponencialmente con la longitud de la secuencia.

Cuando el producto se encoge, que es el caso común con celdas tanh, los primeros pasos de tiempo reciben casi nada de gradiente. La red no puede aprender que algo de hace 80 días importa para el pronóstico, porque la señal de entrenamiento nunca llega tan lejos en el pasado. Este es el problema del gradiente desvanecido, y es la razón por la que la memoria recurrente simple se agota alrededor de unas pocas decenas de pasos. Cuando el producto crece, en cambio, el entrenamiento se vuelve inestable; el recorte de gradientes (gradient clipping) mantiene ese caso bajo control, pero no hay un arreglo igual de simple para el desvanecimiento.

Dos familias de soluciones cambiaron el campo. La LSTM (1997) agrega una celda de memoria con compuertas cuyo estado se actualiza por adición en lugar de multiplicación repetida, dándoles a los gradientes un camino que no se encoge. La atención, introducida por Bahdanau et al. (2014) como un complemento de los modelos recurrentes secuencia-a-secuencia, permite que cada paso de tiempo se conecte directamente con todos los demás; el transformer (Vaswani et al., 2017, «Attention is all you need») eliminó luego la recurrencia por completo, de modo que ningún gradiente tiene que sobrevivir a un producto de 90 pasos. Construimos ambas a continuación.

7. LSTM¶

La celda de memoria de corto-largo plazo (LSTM, Long Short-Term Memory) lleva un estado de celda interno junto al estado oculto y lo controla con tres compuertas aprendidas. La compuerta de olvido decide cuánto del estado de celda anterior conservar, la compuerta de entrada decide cuánta información candidata nueva escribir en él, y la compuerta de salida decide cuánto del estado de celda exponer como estado oculto. Como el estado de celda se actualiza por adición, los gradientes pueden cruzar muchos pasos de tiempo sin desvanecerse. El precio es alrededor de cuatro veces los parámetros de una celda RNN simple del mismo ancho.

Cambiar nn.RNN por nn.LSTM es un cambio de una línea; el resto del modelo y la llamada de entrenamiento son idénticos.

class LSTMForecaster(torch.nn.Module):
    def __init__(self, hidden_size=32, horizon=HORIZON):
        super().__init__()
        self.lstm = torch.nn.LSTM(1, hidden_size, batch_first=True)
        self.head = torch.nn.Linear(hidden_size, horizon)

    def forward(self, x):
        out, _ = self.lstm(x)
        return self.head(out[:, -1, :])


train_model(LSTMForecaster(), "LSTM")

[LSTM] epoch  5  train MSE 0.0047  val MSE 0.0028

[LSTM] epoch 10  train MSE 0.0042  val MSE 0.0025

[LSTM] epoch 15  train MSE 0.0042  val MSE 0.0025

[LSTM] epoch 20  train MSE 0.0042  val MSE 0.0026
[LSTM] 5,470 parameters, trained in 13.6 s, val MAE 1.86 mm

8. Autoatención desde cero¶

La atención toma otra ruta: abandonar la recurrencia y dejar que cada paso de tiempo mire directamente a todos los demás. Cada posición de la ventana emite una consulta (query, «¿qué estoy buscando?»), una clave (key, «¿qué contengo?») y un valor (value, «¿qué transmito?»). La salida en cada posición es un promedio ponderado de todos los valores, con pesos fijados por qué tan bien la consulta de esa posición coincide con cada clave. Una cabeza de atención cabe en unas 15 líneas de PyTorch.

class TinyAttention(torch.nn.Module):
    def __init__(self, d_model=32, horizon=HORIZON):
        super().__init__()
        self.embed = torch.nn.Linear(1, d_model)  # lift each scalar to a vector
        self.Wq = torch.nn.Linear(d_model, d_model, bias=False)
        self.Wk = torch.nn.Linear(d_model, d_model, bias=False)
        self.Wv = torch.nn.Linear(d_model, d_model, bias=False)
        self.head = torch.nn.Linear(d_model, horizon)
        self.scale = d_model ** 0.5

    def forward(self, x):                              # x: (batch, window, 1)
        h = self.embed(x)                              # (batch, window, d_model)
        Q, K, V = self.Wq(h), self.Wk(h), self.Wv(h)   # three views of h
        scores = Q @ K.transpose(1, 2) / self.scale    # (batch, window, window)
        weights = torch.softmax(scores, dim=-1)        # rows sum to one
        context = weights @ V                          # weighted sum of values
        return self.head(context.mean(dim=1))          # mean-pool, then forecast

Línea por línea:

self.embed eleva cada desplazamiento escalar a un vector de dimensión d_model; la atención opera sobre vectores, no sobre escalares.
Wq, Wk, Wv son tres mapas lineales que producen las consultas Q, las claves K y los valores V, cada uno de forma (batch, 90, d_model). Son tres «vistas» aprendidas de la misma secuencia embebida.
Q @ K.transpose(1, 2) calcula todos los productos punto consulta-clave de una vez: una matriz de puntajes de 90 por 90 por muestra. La entrada $(i, j)$ mide qué tan relevante es el día $j$ para el día $i$ .
Dividir por $\sqrt{d_{model}}$ mantiene los puntajes cerca de la escala unitaria para que la softmax no se sature.
softmax(scores, dim=-1) convierte cada fila de puntajes en pesos positivos que suman uno.
weights @ V forma el promedio ponderado: la salida de cada día mezcla información de los 90 días, cercanos o lejanos, al mismo costo. Ningún producto de jacobianas de 90 pasos en ninguna parte.
Promediamos sobre las 90 posiciones (mean pooling) y aplicamos una cabeza lineal para producir los 30 valores de pronóstico.

Note lo que falta: nada en estas líneas conoce el orden de los pasos de tiempo. Baraje los 90 días de una ventana y la salida promediada es idéntica. El transformer corrige esto con codificaciones posicionales.

train_model(TinyAttention(), "Attention")

[Attention] epoch  5  train MSE 0.0047  val MSE 0.0028

[Attention] epoch 10  train MSE 0.0045  val MSE 0.0026

[Attention] epoch 15  train MSE 0.0045  val MSE 0.0026

[Attention] epoch 20  train MSE 0.0045  val MSE 0.0026
[Attention] 4,126 parameters, trained in 8.5 s, val MAE 1.87 mm

9. Codificador transformer¶

Una capa de codificador transformer es autoatención más una pequeña red feed-forward, con conexiones residuales y normalización de capa alrededor de cada una, y normalmente varias cabezas de atención en paralelo. PyTorch trae el bloque completo como torch.nn.TransformerEncoderLayer; apilamos dos de ellos con torch.nn.TransformerEncoder.

Como la atención es ciega al orden, primero agregamos una codificación posicional a las entradas embebidas: un patrón fijo de senos y cosenos a distintas frecuencias, un vector por posición. Después de esta suma, el día 3 y el día 73 de una ventana se ven distintos para el modelo aun cuando sus valores de desplazamiento sean iguales.

class PositionalEncoding(torch.nn.Module):
    def __init__(self, d_model, max_len=500):
        super().__init__()
        pos = torch.arange(max_len).unsqueeze(1).float()
        freq = 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 * freq)
        pe[:, 1::2] = torch.cos(pos * freq)
        self.register_buffer("pe", pe)

    def forward(self, x):                 # x: (batch, seq, d_model)
        return x + self.pe[: x.shape[1]]


class TinyTransformer(torch.nn.Module):
    def __init__(self, d_model=32, nhead=4, num_layers=2, horizon=HORIZON):
        super().__init__()
        self.embed = torch.nn.Linear(1, d_model)
        self.pos = PositionalEncoding(d_model)
        layer = torch.nn.TransformerEncoderLayer(
            d_model=d_model, nhead=nhead, dim_feedforward=64,
            dropout=0.0, batch_first=True)
        self.encoder = torch.nn.TransformerEncoder(
            layer, num_layers=num_layers, enable_nested_tensor=False)
        self.head = torch.nn.Linear(d_model, horizon)

    def forward(self, x):
        h = self.encoder(self.pos(self.embed(x)))
        return self.head(h.mean(dim=1))


train_model(TinyTransformer(), "Transformer")

[Transformer] epoch  5  train MSE 0.0049  val MSE 0.0031

[Transformer] epoch 10  train MSE 0.0049  val MSE 0.0031

[Transformer] epoch 15  train MSE 0.0046  val MSE 0.0027

[Transformer] epoch 20  train MSE 0.0042  val MSE 0.0025
[Transformer] 18,142 parameters, trained in 23.8 s, val MAE 1.84 mm

10. Comparación¶

Los cuatro modelos se entrenaron con las mismas ventanas y se calificaron con la misma métrica, así que los números de abajo son directamente comparables. La tabla también lleva dos filas que no cuestan nada entrenar:

Persistencia: cada día del pronóstico es igual al último día del contexto. En nuestras coordenadas ancladas ese pronóstico es exactamente cero, así que su MAE es una línea de código.
Ingenuo estacional: cada día pronosticado es igual al valor observado 365 días antes — el modelo de referencia correcto para una serie puramente estacional.

Un pronosticador aprendido que no vence a la persistencia no ha aprendido nada; esa es la regla de lectura para cada fila de arquitectura de abajo.

# Baselines that require no training
# persistence: in anchored coordinates the forecast is identically zero
mae_persistence = Y_val.abs().mean().item() * sigma

# seasonal naive: each forecast day equals the value 365 days earlier
starts = np.arange(split, n_days - WINDOW - HORIZON + 1)  # anchor day of each validation pair
seas_pred = np.stack([z[i + WINDOW - 365 : i + WINDOW + HORIZON - 365] for i in starts])
true_z = Y_val.numpy() + anc_val.numpy()                  # targets back in absolute z units
mae_seasonal = float(np.abs(seas_pred - true_z).mean() * sigma)

baseline_rows = [
    {"model": "persistence (last value)", "parameters": 0, "train time (s)": 0.0,
     "val MAE, 30-day horizon (mm)": round(mae_persistence, 2)},
    {"model": "seasonal naive (365 d earlier)", "parameters": 0, "train time (s)": 0.0,
     "val MAE, 30-day horizon (mm)": round(mae_seasonal, 2)},
]
comparison = pd.DataFrame(
    baseline_rows +
    [{"model": name,
      "parameters": r["n_params"],
      "train time (s)": round(r["train_time_s"], 1),
      "val MAE, 30-day horizon (mm)": round(r["val_mae_mm"], 2)}
     for name, r in results.items()]
).set_index("model")
comparison

model_colors = {"RNN": "#4477AA", "LSTM": "#EE6677",
                "Attention": "#228833", "Transformer": "#AA3377"}

fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))

for name, r in results.items():
    axes[0].plot(r["history"]["val"], color=model_colors[name], label=name)
axes[0].set_yscale("log")
axes[0].set_xlabel("epoch")
axes[0].set_ylabel("validation MSE (standardized units)")
axes[0].set_title("Validation loss")
axes[0].legend()
axes[0].grid(alpha=0.3)

i = 300  # same validation window as before
a = anc_val[i].item()
t_ctx = np.arange(WINDOW)
t_fut = np.arange(WINDOW, WINDOW + HORIZON)
axes[1].plot(t_ctx, (X_val[i, :, 0].numpy() + a) * sigma + mu,
             color="#333333", lw=1, label="context")
axes[1].plot(t_fut, (Y_val[i].numpy() + a) * sigma + mu,
             color="#333333", lw=2, ls="--", label="truth")
for name, r in results.items():
    r["model"].eval()
    with torch.no_grad():
        pred = r["model"](X_val[i : i + 1].to(device)).cpu().numpy()[0]
    axes[1].plot(t_fut, (pred + a) * sigma + mu, color=model_colors[name],
                 lw=1.2, label=name)
axes[1].set_xlabel("day within window")
axes[1].set_ylabel("displacement (mm)")
axes[1].set_title("One validation forecast")
axes[1].legend(fontsize=8)
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

Empiece por las filas de referencia. La persistencia marca 2.10 mm, y los cuatro modelos aprendidos quedan en 1.83-1.87 mm: una mejora de alrededor del 13 %, real pero modesta, que es el titular correcto para una serie tan dominada por una tendencia suave y la estacionalidad. El ingenuo estacional fracasa con 15.5 mm — el valor de 365 días antes queda sistemáticamente unos 12 mm por debajo con una tendencia tectónica de 12 mm/año —, un recordatorio de que un modelo de referencia solo informa cuando coincide con la estructura a la que apunta. Los cuatro puntajes aprendidos quedan muy cerca entre sí, y eso es lo esperable. Esta serie sintética está dominada por una tendencia y un ciclo estacional que cualquier arquitectura puede extraer de un contexto de 90 días, y 20 épocas sobre unos pocos miles de ventanas dejan ruido en el ordenamiento; vuelva a correr con otra semilla y el ranking puede barajarse. El ranking no es la lección. La mecánica sí lo es: cómo mueve cada arquitectura la información a través del tiempo, y qué cuesta eso en parámetros y en comportamiento del gradiente.

Ejercicio¶

Duplique el horizonte de pronóstico a 60 días (fije HORIZON = 60, reconstruya las ventanas, reentrene los cuatro modelos). ¿Cómo cambian los MAE de validación, y por qué?

Solución

Los cuatro MAE crecen, porque la incertidumbre se acumula con la anticipación: el día 60 depende del ruido, de la fase estacional y de cualquier comportamiento transitorio que el contexto de 90 días restringe cada vez menos. El error tampoco es uniforme a lo largo del horizonte; si calcula el MAE por día de anticipación, los primeros días del pronóstico son mucho mejores que los últimos. El ranking relativo de los cuatro modelos suele quedar dentro del ruido en esta serie sintética.

11. Resumen¶

Una tarea de pronóstico se define por su ventana de contexto y su horizonte de pronóstico; los pares supervisados salen de deslizar esa ventana sobre la serie, con una división temporal entrenamiento/validación para que ninguna información se filtre a través de la frontera.
Una RNN simple lleva un estado oculto a través del tiempo pero no puede aprender estructura de largo alcance, porque la BPTT multiplica aproximadamente una jacobiana por paso de tiempo y el producto se desvanece o explota.
La LSTM encamina la información a través de un estado de celda actualizado por adición, controlado por las compuertas de olvido, entrada y salida, de modo que los gradientes sobreviven a secuencias largas.
La autoatención conecta cada paso de tiempo con todos los demás en una sola operación: consultas, claves, valores, productos punto escalados, softmax, suma ponderada. Un codificador transformer envuelve eso en conexiones residuales, normalización de capa y bloques feed-forward, más codificaciones posicionales para restaurar la información de orden.

La atención es el bloque constructivo detrás de los sistemas de pronóstico actuales, desde los modelos meteorológicos de aprendizaje automático hasta los pronosticadores de series de tiempo preentrenados. El cuaderno 4.10 corre una comparación de pronóstico de estas arquitecturas sobre datos geocientíficos reales.