XII.2 - Travaux pratiques - Prédiction de séries temporelles

L’objectif de cette séance de travaux pratiques est de mettre en œuvre un réseau de neurones récurrent pour prédire les variations d’une série temporelle. En particulier, nous entraînerons un modèle de régression qui pourra servir soit à régresser des données manquantes, soit à prévoir les valeurs futures de la série. L’illustration se fera sur un jeu de données météorologiques.

Création d’un jeu d’entraînement et de validation

Nous allons nous entraîner sur un jeu de données météorologiques disponible à l’adresse suivante : https://cedric.cnam.fr/vertigo/Cours/ml2/docs/meteo_dataset.csv.

Pour commencer, téléchargeons ce jeu de données.

import os
import urllib.request

if not os.path.exists('meteo_dataset.csv'):
    urllib.request.urlretrieve('https://cedric.cnam.fr/vertigo/Cours/ml2/docs/meteo_dataset.csv', 'meteo_dataset.csv')

Il s’agit d’un jeu de données au format CSV, que nous pouvons visualiser, par exemple, avec la bibliothèque pandas (https://pandas.pydata.org/) en utilisant la méthode pandas.read_csv.

import pandas as pd

data = pd.read_csv('meteo_dataset.csv')
data

Comme nous pouvons le voir sur la sortie ci-dessus, le jeu de données meteo est une série de \(N = 420551\) observations. Chaque observation est décrite par un ensemble de variables, des grandeurs physiques issues d’un système météorologique. Formellement, nous travaillons sur une séquence de vecteurs \((\mathbf u_n)_{0\leq n \leq N-1} \in \mathbb R^{7}\). Les différentes composantes de ces vecteurs correspondent aux paramètres physiques mesurés : \(\mathcal P = \{\)pression, température, saturation, déficit de pression de vapeur, humidité, densité, vitesse du vent\(\}\).

Nous pouvons observer les variations de certains de ces paramètres au cours du temps pour l’ensemble du jeu de données, en traçant leur évolution au cours du temps avec matplotlib :

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(20,7))
# Évolution de la pression
ax1 = fig.add_subplot(221)
plt.plot(data['Pressure'])
ax1.title.set_text('Pression')
# Évolution de la température
ax2 = fig.add_subplot(222)
plt.plot(data['Temperature'])
ax2.title.set_text('Température')
# Évolution de la saturation
ax3 = fig.add_subplot(223)
plt.plot(data['Saturation vapor pressure'])
ax3.title.set_text('Saturation')
# Évolution du déficit de pression de vapeur
ax4 = fig.add_subplot(224)
plt.plot(data['Vapor pressure deficit'])
ax4.title.set_text('Déficit de pression de vapeur')

plt.show()

Question

Commenter la forme des courbes. Que peut-on dire d’après l’aspect global des données ? Parmi les quatre paramètres que nous avons affiché, lesquels sont susceptibles d’être plus difficiles à apprendre que les autres ?

Correction

La température, la saturation et le déficit de pression de vapeur ont des évolutions périodiques très marquées, avec une cohérence temporelle élevée. La courbe de température est notamment assez peu bruitée. À l’inverse, la courbe de pression est bruitée, avec des outliers importants et ne semble pas exhiber de périodicité évidente. Ce paramètre risque donc d’être nettement plus difficile à prédire.

Division du jeu de données en apprentissage et validation

Comme d’habitude, nous allons séparer les données en un ensemble d’apprentissage et un ensemble de validation. L’ensemble de validation nous servira à l’évaluation du modèle. Dans notre cas, nous choisissons arbitrairement d’utiliser 80% des données pour l’apprentissage et les 20% restants pour la validation.

Question

Quelles précautions doit-on prendre dans la division de ce jeu de données en ensemble de train et de test ? Pourquoi ?

Correction

Comme il s’agit de séries temporelles, les observations ne sont pas i.i.d. (indépendantes et identiquement distribuées). En effet, l’observation \(x_t\) au pas de temps \(t\) sera très certainement corrélée avec l’observation \(x_{t+1}\) au pas de temps suivant. Il ne faut donc pas réaliser une division aléatoire avec un tirage uniforme : les données de test seraient très similaires aux données d’entraînement et notre estimation de l’erreur de généralisation sera bien trop optimiste.

Une solution consiste à entraîner sur les premières valeurs de la série temporelle, puis de s’évaluer sur les dernières valeurs, en créant un jeu d’apprentissage et un jeu de test qui sont temporellement disjoints.

split_fraction = 0.80
split_idx = int(split_fraction * len(data))

# Train: on conserve les 80% premières lignes
train_data = data.loc[0:split_idx-1].values
# Test: on conserve les 20% dernières lignes
val_data = data.loc[split_idx:].values

print(f"Jeu d'apprentissage : {train_data.shape}, jeu d'évaluation : {val_data.shape}")

Dans un premier temps, nous allons nous intéresser à la prédiction de la température à partir des six autres paramètres de la série temporelle. Il s’agit donc d’une configuration many-to-many, c’est-à-dire que le modèle décisionnel consiste à prédire une séquence \(\mathbf y_t\) (ici, la température) à partir d’une séquence de variables explicatives \(\mathbf x_t\) (les autres paramètres).

Question

Dans cette configuration, quelles sont les données ? Quelles sont les étiquettes ? Quelles sont leurs dimensions respectives ?

Correction

Le modèle prendra en entrée une séquence de \(M\) éléments correspondant à tous les paramètres sauf la température pour \(M\) acquisitions et prédira pour chaque pas de temps la température associée. Formellement, le modèle prendra en entrée une séquence \(\mathbf X_n = [\mathbf x_{n},\cdots,\mathbf x_{n+M-1}] \in \mathbb R^{M\times 6}\) telle que \(\forall\, n \in \{0,\cdots , N-1\}, p \in \mathcal P\setminus \{\text{Temperature}\}, \mathbf x_{n,p} = \mathbf u_{n,p}\). En sortie, le modèle prédit la séquence \(\mathbf Y_n = [ y_{n},\cdots, y_{n+M-1}] \in \mathbb R^{M\times 1}\) telle que \(\forall\, n \in \{0,\cdots , N-1\}, y_{n} = \mathbf u_{n,\text{Temperature}}\).

Afin d’entraîner notre réseau de neurones par descente de gradient, nous allons devoir subdiviser le jeu de données en batch de séquences. En PyTorch, il n’existe pas d’équivalent direct de la fonction timeseries_dataset_from_array de Keras : nous allons donc écrire nous-mêmes un petit Dataset qui, à partir du tableau de la série temporelle, découpe des fenêtres glissantes de longueur sequence_length.

import torch
from torch.utils.data import Dataset, DataLoader

# Graine aléatoire, pour la reproductibilité (ordre des batches, initialisation du GRU)
torch.manual_seed(42)

class TimeseriesDataset(Dataset):
    """ Convertit un tableau (N, d) en fenêtres glissantes de longueur `sequence_length`,
        échantillonnées tous les `sampling_rate` pas de temps, pour prédire à chaque pas de
        temps de la fenêtre la valeur correspondante de `targets` (many-to-many). """
    def __init__(self, data, targets, sequence_length, sampling_rate):
        self.data = torch.as_tensor(data, dtype=torch.float32)
        self.targets = torch.as_tensor(targets, dtype=torch.float32)
        self.sequence_length = sequence_length
        self.sampling_rate = sampling_rate
        self.window_span = sequence_length * sampling_rate

    def __len__(self):
        return len(self.data) - self.window_span + 1

    def __getitem__(self, idx):
        x = self.data[idx: idx + self.window_span: self.sampling_rate]
        y = self.targets[idx: idx + self.window_span: self.sampling_rate]
        return x, y

Question

En vous aidant du code ci-dessus, expliquer à quoi correspondent les paramètres sequence_length et sampling_rate.

Correction

sequence_length correspond à la taille de la fenêtre de la série temporelle, c’est-à-dire la longueur maximale de la séquence considérée. sampling_rate définit le taux d’échantillonnage et permet de réduire la dimension de la série en considérant par exemple qu’une valeur sur deux.

Question

Compléter le code suivant pour créer un jeu de données d’entraînement au format attendu par PyTorch (à l’aide de la classe TimeseriesDataset définie ci-dessus). Répéter pour le jeu de données d’évaluation. On se fixera des batch de taille \(256\), des séquences de longueur \(120\) éléments et, pour réduire le temps de calcul, nous allons sous-échantillonner les données d’un facteur \(6\).

sampling_rate = 6 # Fréquence d'échantillonnage
sequence_length = 120  # Longueur de la séquence
batch_size = 256 # Taille de batch

# Création du jeu d'entraînement

dataset_train = TimeseriesDataset(
    data=train_data[:, (0,2,3,4,5,6)], # sélectionne toutes les colonnes sauf la température
    targets=train_data[:, 1],           # sélectionne la colonne de température
    # À compléter
)
train_loader = DataLoader(dataset_train, batch_size=batch_size, shuffle=True)

# Création du jeu d'évaluation

dataset_val = ... # À compléter
val_loader = ... # À compléter

Correction

sampling_rate = 6 # Fréquence d'échantillonnage
sequence_length = 120  # Longueur de la séquence
batch_size = 256 # Taille de batch

# Création du jeu d'entraînement

dataset_train = TimeseriesDataset(
    data=train_data[:, (0,2,3,4,5,6)],
    targets=train_data[:, 1],
    sequence_length=sequence_length,
    sampling_rate=sampling_rate,
)
train_loader = DataLoader(dataset_train, batch_size=batch_size, shuffle=True)

# Création du jeu d'évaluation

dataset_val = TimeseriesDataset(
    data=val_data[:, (0,2,3,4,5,6)],
    targets=val_data[:, 1],
    sequence_length=sequence_length,
    sampling_rate=sampling_rate,
)
val_loader = DataLoader(dataset_val, batch_size=batch_size, shuffle=False)

Entraînement du réseau récurrent

Maintenant que le jeu de données est mis en forme, nous pouvons définir le réseau récurrent que nous allons utiliser pour prédire la température à partir des autres paramètres météorologiques. Pour ce TP, nous allons implémenter un réseau récurrent simple constitué d’une Gated Recurrent Unit (GRU) et de deux couches linéaires.

from torch import nn

Question

Quelles sont les dimensions des vecteurs en entrée et en sortie du réseau ?

Question

Compléter le code ci-dessous qui réalise la construction du modèle, sous la forme d’une classe héritant de nn.Module. On choisira d’utiliser une dimension égale à celle des variables d’entrée pour la couche cachée du GRU. On rappelle que nn.GRU renvoie, pour une séquence d’entrée, la séquence complète des états cachés (la sortie du GRU à chaque pas de temps), de dimensions (batch, sequence_length, hidden_size).

class GRUModel(nn.Module):
    def __init__(self, input_dim):
        super().__init__()
        # À compléter
        ...

    def forward(self, x):
        # À compléter
        ...

dimension = train_data.shape[1] - 1
model = GRUModel(dimension)

Correction

class GRUModel(nn.Module):
    def __init__(self, input_dim):
        super().__init__()
        self.gru = nn.GRU(input_dim, input_dim, batch_first=True)
        self.dense = nn.Linear(input_dim, input_dim)
        self.relu = nn.ReLU()
        self.output = nn.Linear(input_dim, 1)

    def forward(self, x):
        gru_out, _ = self.gru(x)
        dense_out = self.relu(self.dense(gru_out))
        return self.output(dense_out)

dimension = train_data.shape[1] - 1
model = GRUModel(dimension)

Note

batch_first=True indique à nn.GRU que la dimension du batch vient en premier dans les tenseurs d’entrée et de sortie (forme (batch, séquence, features)), comme c’est le cas pour la plupart des autres couches PyTorch. Par défaut, nn.GRU attend au contraire la séquence en premier ((séquence, batch, features)), pour des raisons historiques de performance.

Concernant l’optimisation, nous allons utiliser le solveur Adam pour la descente de gradient. On fixe le pas d’apprentissage à \(0,001\) pour l’instant. On utilisera l’erreur quadratique moyenne (nn.MSELoss) comme fonction de coût, adaptée à un problème de régression :

learning_rate = 0.001
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

Comme pour les séances précédentes, nous devons cette fois écrire nous-mêmes la boucle d’apprentissage. Nous allons entraîner ce modèle pendant 5 époques, et suivre au passage l’erreur de validation :

epochs = 5

for epoch in range(epochs):
    model.train()
    for x, y in train_loader:
        y_pred = model(x).squeeze(-1)
        loss = criterion(y_pred, y)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

    model.eval()
    with torch.no_grad():
        val_losses = [criterion(model(x).squeeze(-1), y).item() for x, y in val_loader]
    print(f"Epoch {epoch}: loss (train) = {loss.item():.4f}, loss (val) = {sum(val_losses)/len(val_losses):.4f}")

Note

model.train() et model.eval() jouent, pour un modèle PyTorch personnalisé, le même rôle que nous avons déjà vu : ils permettent aux couches sensibles au mode d’exécution (dropout, batch normalization…) de se comporter différemment en apprentissage et en évaluation. Notre GRU n’en utilise pas, mais c’est une bonne habitude à prendre systématiquement.

L’affichage de l’erreur de validation à chaque époque nous permet de vérifier que nous ne sommes pas en régime de sur-apprentissage. Pour une analyse plus qualitative, nous pouvons visualiser à présent les résultats des prédictions sur les 5 premiers éléments du jeu de validation. La fonction ci-dessous show_plot permet de visualiser sur un même graphique la prédiction et la vérité terrain.

def show_plot(predicted_data, true_data, title=None):
    plt.title(title)
    plt.plot(predicted_data, "rx", label="Température prédite")
    plt.plot(true_data, ".-", label="Vraie température")
    plt.legend()
    plt.xlabel("Pas de temps")
    plt.show()
    return

Question

À l’aide du modèle entraîné, calculer les prédictions de température pour les trois premières séquences x du jeu de données de validation. Utiliser la fonction show_plot ci-dessus pour afficher les graphiques correspondants. On pensera à repasser le modèle en mode évaluation (model.eval()) et à désactiver le calcul des gradients.

model.eval()
with torch.no_grad():
    for x, y in list(val_loader)[:1]:
        # À compléter, sur les 3 premières séquences du batch
        ...

Correction

model.eval()
with torch.no_grad():
    x, y = next(iter(val_loader))
    y_pred = model(x).squeeze(-1)
    for i in range(3):
        show_plot(y_pred[i].numpy(), y[i].numpy(), title="Prédiction synchrone")

Question

Que constatez-vous concernant les performances du modèle ?

Correction

Les prédictions de température sont très proches de la vérité terrain. Les performances en train et en test sont bonnes et à peu près équivalentes, ce qui signifie a priori qu’il n’y a pas de surapprentrissage.

Bien que simple, notre modèle récurrent semble à même de prédire avec une haute précision la température à partir des autres paramètres météorologiques. Même si nous avons utilisé l’erreur quadratique moyenne pour l’optimisation du modèle, nous pouvons calculer l’erreur absolue moyenne, qui est souvent plus facile à interpréter :

mean_absolute_error = nn.L1Loss()
error = 0

model.eval()
with torch.no_grad():
    for x, y in val_loader:
        y_pred = model(x).squeeze(-1)
        error += mean_absolute_error(y_pred, y).item()

print(f"Erreur absolue moyenne : {error/len(val_loader):.5f}")

Question

À faire chez vous pour approfondir, si vous avez le temps. Refaire l’expérience mais en prédisant la pression à partir des autres paramètres. Les conclusions sont-elles identiques ?

Prédiction de la température future

Nous avons vu précédemment comment prédire une séquence de valeurs à partir d’une autre séquence. Pour aller plus loin, nous allons à présent chercher non plus à reconstruire température à partir des autres paramètres, mais à prédire la température future, c’est-à-dire à \(t+1\), à partir de toutes les variables météorologiques jusqu’au temps \(t\).

Comme le but est de prédire la température future, nous allons diviser nos données en séquences de \(M\) éléments. Chacune aura pour étiquette associée la température \(T\) instants plus tard. Formellement, notre modèle prendra en entrée une séquence décrite par \(\mathbf X_n = \{\mathbf u_{n}, \cdots, \mathbf u_{n+M}\} \in \mathbb R^{M\times 7}\) et prédira en sortie un scalaire unique \(y_n = \mathbf u_{n+M+T+1,\text{Temperature}} \in \mathbb R\). Il s’agit donc désormais d’un problème many-to-one.

Nous utiliserons une séquence de \(M = 720\) observations passées pour prédire la température à \(T = 6\) instants plus tard. Les séquences d’entrée démarreront donc à l’indice \(0\) tandis que la première étiquette correspondra à l’indice \(720+6=726\). Ici encore, nous allons sous-échantillonner la séquence d’origine par un facteur \(6\), ce qui nous donnera en entrée des séquences de longueur \(720/6=120\), comme précédemment. Avec ce sous-échantillonnage, cela revient donc à faire la prédiction \(t+1\) à partir des \(120\) valeurs précédentes.

Cette fois, chaque séquence d’entrée est associée à une unique étiquette (et non plus une séquence d’étiquettes) : nous avons besoin d’un second petit Dataset, adapté à ce cas many-to-one.

class FutureTimeseriesDataset(Dataset):
    """ Fenêtres glissantes de longueur `sequence_length` (échantillonnées tous les
        `sampling_rate` pas de temps), chacune associée à une unique étiquette (celle de
        `targets` à `future_step` pas de temps après la fin de la fenêtre). """
    def __init__(self, data, targets, sequence_length, sampling_rate, future_step):
        self.data = torch.as_tensor(data, dtype=torch.float32)
        self.targets = torch.as_tensor(targets, dtype=torch.float32)
        self.sequence_length = sequence_length
        self.sampling_rate = sampling_rate
        self.window_span = sequence_length * sampling_rate
        self.future_step = future_step

    def __len__(self):
        return len(self.data) - self.window_span - self.future_step

    def __getitem__(self, idx):
        x = self.data[idx: idx + self.window_span: self.sampling_rate]
        y = self.targets[idx + self.window_span + self.future_step]
        return x, y

past_length = 720
future_step = 6
step_size = 6 # Facteur d'échantillonnage
sequence_length = int(past_length / step_size)
batch_size = 256

dataset_train = FutureTimeseriesDataset(
    data=train_data,
    targets=train_data[:, 1],
    sequence_length=sequence_length,
    sampling_rate=step_size,
    future_step=future_step,
)
train_loader = DataLoader(dataset_train, batch_size=batch_size, shuffle=True)

dataset_val = FutureTimeseriesDataset(
    data=val_data,
    targets=val_data[:, 1],
    sequence_length=sequence_length,
    sampling_rate=step_size,
    future_step=future_step,
)
val_loader = DataLoader(dataset_val, batch_size=batch_size, shuffle=False)

Correction

Comme précédemment, nous allons entraîner un réseau récurrent constitué d’une couche récurrente utilisant la cellule GRU et de deux couches linéaires séparées par une activation ReLU.

Question

Quelle est la taille des éléments en entrée et en sortie du réseau ?

Question

Compléter le code ci-dessous qui gère la définition du modèle. Cette fois, seul le dernier état caché du GRU nous intéresse (puisque l’on prédit un scalaire unique, pas une séquence) : on pourra utiliser le second élément renvoyé par self.gru(x), de dimensions (1, batch, hidden_size) (à ne pas confondre avec la séquence complète des états cachés). On utilisera une couche cachée de dimension 32. Ne pas oublier de créer un optimiseur Adam avec un pas d’apprentissage de \(0,001\). On utilisera l’erreur quadratique moyenne (MSE) comme fonction de coût.

class GRUFutureModel(nn.Module):
    def __init__(self, input_dim, hidden_dim):
        super().__init__()
        # À compléter
        ...

    def forward(self, x):
        # À compléter
        ...

input_dim = train_data.shape[1]
model = GRUFutureModel(input_dim, hidden_dim=32)
criterion = ... # À compléter
optimizer = ... # À compléter

Correction

class GRUFutureModel(nn.Module):
    def __init__(self, input_dim, hidden_dim):
        super().__init__()
        self.gru = nn.GRU(input_dim, hidden_dim, batch_first=True)
        self.dense = nn.Linear(hidden_dim, hidden_dim)
        self.relu = nn.ReLU()
        self.output = nn.Linear(hidden_dim, 1)

    def forward(self, x):
        _, h_n = self.gru(x)   # h_n : dernier état caché, (1, batch, hidden_dim)
        h_n = h_n.squeeze(0)   # (batch, hidden_dim)
        dense_out = self.relu(self.dense(h_n))
        return self.output(dense_out).squeeze(-1)

learning_rate = 0.001
input_dim = train_data.shape[1]
hidden_dim = 32
model = GRUFutureModel(input_dim, hidden_dim)
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

À nouveau, nous pouvons entraîner le modèle sur le jeu de données à l’aide de la boucle d’apprentissage habituelle :

epochs = 5

for epoch in range(epochs):
    model.train()
    for x, y in train_loader:
        y_pred = model(x)
        loss = criterion(y_pred, y)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

    model.eval()
    with torch.no_grad():
        val_losses = [criterion(model(x), y).item() for x, y in val_loader]
    print(f"Epoch {epoch}: loss (train) = {loss.item():.4f}, loss (val) = {sum(val_losses)/len(val_losses):.4f}")

Pour terminer, nous pouvons désormais visualiser les résultats en prédiction de notre modèle sur les trois premières séquences du jeu de validation :

def show_plot(sequence, prediction, ground_truth, delta=future_step/step_size, title=None):
    plt.title(title)
    plt.plot(sequence.flatten(), ".-", label="Température passée")
    plt.plot(len(sequence) + delta, prediction, "go", label="Prédiction")
    plt.plot(len(sequence) + delta, ground_truth, "rx", label="Vraie température")
    plt.legend()
    plt.xlim(-1, len(sequence) + delta + 3)
    plt.xlabel("Pas de temps")
    plt.show()

model.eval()
with torch.no_grad():
    x, y = next(iter(val_loader))
    y_pred = model(x)
    for i in range(5):
        show_plot(x[i][:, 1].numpy(), y_pred[i].item(), y[i].item(), title="Prédiction de la température")

Question

Qu’observez-vous qualitativement dans les prédictions de notre modèle ? Que pouvez-vous en conclure sur la difficulté de la tâche de prédiction par rapport à la tâche de régression précédemment considérée ?

Correction

Les prédictions sont moins souvent correctes que précédemment. Il y a plus de grandes erreurs de prédictions entre la température prédite et la température future réelle. Ce n’est pas inattendu car l’on pouvait s’attendre à ce que la prédiction des valeurs futures soit un problème plus difficile que la reconstruction de la température à un pas de temps donné en connaissant toutes les autres variables.