2.6 Méthodes de rééchantillonnage - Machine Learning en géosciences

Les données ont une variabilité naturelle. Toute statistique calculée sur un échantillon (une moyenne, une corrélation, une pente de régression) hérite de cette variabilité.

Les méthodes de rééchantillonnage permettent de la mesurer. Au sens large, le rééchantillonnage désigne toute technique où l’on tire de façon répétée des observations d’un échantillon, où l’on recalcule une statistique sur chaque tirage et où l’on étudie la distribution des résultats. Les applications incluent les tests d’hypothèse, la propagation d’incertitude et les intervalles de confiance.

La leçon comporte trois niveaux. Le niveau 1 parcourt trois techniques de rééchantillonnage de base : la randomisation, le bootstrap et Monte-Carlo. Le niveau 2 applique le bootstrap à l’inférence de modèle pour une régression linéaire, d’abord sur des données GNSS synthétiques dont la réponse est connue, puis sur des données GNSS réelles de la zone de subduction de Cascadia. Le niveau 3 se tourne vers l’autre sens du mot rééchantillonnage — changer l’échantillonnage d’un signal : sous-échantillonner sans repliement (aliasing), interpoler les lacunes selon une politique explicite, agréger des réseaux de stations irréguliers, et le bootstrap par blocs pour le bruit corrélé.

D’abord, importons les modules nécessaires :

🖥️ Diapositives du cours — Séance 06 (lun. 12 oct.)

import os
import requests
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import mlgeo_synth

Une note sur les nombres aléatoires. Le code NumPy ancien initialise l’état aléatoire global avec np.random.seed(42) puis appelle des fonctions comme np.random.normal. Le NumPy moderne remplace cela par un objet générateur explicite : rng = np.random.default_rng(42). Le générateur porte son propre état, si bien que différentes parties d’un programme (ou d’un carnet) n’interfèrent pas entre elles à travers un état global caché. Nous utilisons l’idiome du générateur dans tout ce livre. Fixer la graine rend le carnet reproductible de bout en bout.

# One generator for the whole notebook, with a fixed seed for reproducibility.
rng = np.random.default_rng(42)
rng

Generator(PCG64) at 0x14F424820

1. Exemples de techniques de rééchantillonnage (niveau 1)¶

1.1 Randomisation¶

Étant donné deux jeux de données, $A$ et $B$ , et un paramètre $x$ , nous pouvons réaffecter aléatoirement les observations à $A$ ou à $B$ , calculer une statistique (par exemple $x_A - x_B$ ) et répéter pour construire une distribution de cette statistique sous l’hypothèse nulle que les étiquettes de groupe n’ont pas d’importance.

# We begin with two datasets, A and B
A = rng.normal(5, 2.5, 100)
B = rng.normal(5.5, 2.5, 100)

# We then calculate the means of each dataset
mean_A = np.mean(A)
mean_B = np.mean(B)

print(f'The means of A and B are {mean_A:.3f} and {mean_B:.3f}, respectively.')

# And, for the sake of illustration, also calculate the difference between these means
diff_means = mean_A - mean_B
print(f'The difference of means is {diff_means:.3f}.')

The means of A and B are 4.874 and 5.473, respectively.
The difference of means is -0.599.

Maintenant, rééchantillonnons.

# First, how many times do we want to resample?
number_runs = 10000

# Next, we create an array that will store the difference of means
array_of_diffs = np.zeros(number_runs)

# To ease computational burden, declare some variables:
# a combined list of A and B, and the length of A
combined = np.concatenate((A, B))
length_A = len(A)

À chaque tirage, nous mélangeons les données combinées, les scindons en deux nouveaux groupes des tailles d’origine et recalculons la différence des moyennes.

for i in range(number_runs):
    # Shuffle the combined list. Note that shuffle works in place!
    rng.shuffle(combined)
    # Split the list into A and B, maintaining their original sizes.
    new_A = combined[0:length_A]
    new_B = combined[length_A:len(combined)]
    # Calculate and store a difference of means
    array_of_diffs[i] = np.mean(new_A) - np.mean(new_B)

# Plot the array of diffs
plt.hist(array_of_diffs, color='black')

# Given an alpha of 0.05, can we accept or reject the null hypothesis of no difference in means?
alpha = 0.05
lower_critical_value = np.quantile(array_of_diffs, alpha / 2)
upper_critical_value = np.quantile(array_of_diffs, 1 - (alpha / 2))

# Plot the critical values and the observed value
plt.axvline(x=diff_means, color='r')
plt.axvline(x=lower_critical_value, color='g')
plt.axvline(x=upper_critical_value, color='b')
plt.xlabel('Difference of means')
plt.ylabel('Count')
plt.legend(['Observed difference of means', 'Lower critical value', 'Upper critical value'],
           loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()

1.2 Le bootstrap¶

Avec le bootstrap, on tire de façon répétée des observations, avec remise, d’un échantillon. Chaque rééchantillon bootstrap a la même taille que l’échantillon d’origine. Le tirage étant avec remise, un rééchantillon contiendra des observations en double et en omettra d’autres. La dispersion de la statistique à travers les rééchantillons estime l’incertitude d’échantillonnage de cette statistique.

Le bootstrap ne fait pas d’hypothèses fortes sur la distribution sous-jacente des données.

Nous créons d’abord une « vraie population » synthétique de données corrélées avec la méthode multivariate_normal du générateur (documentation ici). La moyenne des deux variables est nulle et leur covariance vaut -0,75.

# A population where two variables are strongly anticorrelated
correlated_data = rng.multivariate_normal([0, 0], [[1, -0.75], [-0.75, 1]], 1000)

Vérifiez que les données sont bien anticorrélées en traçant une variable contre l’autre.

plt.scatter(correlated_data[:, 0], correlated_data[:, 1], marker='x', c='black')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(['Observation'])
plt.show()

Nous vérifions que le coefficient de corrélation de Pearson est proche de notre cible avec la fonction NumPy corrcoef.

correlation_matrix = np.corrcoef(correlated_data[:, 0], correlated_data[:, 1])
print(f'The population correlation coefficient is {correlation_matrix[0, 1]:.3f}.')

The population correlation coefficient is -0.721.

En pratique, nous n’observons jamais la population entière. Nous prenons maintenant un petit sous-ensemble des données — voyez ce sous-ensemble comme l’échantillon que nous avons réellement collecté.

nsubset = 50
subset = rng.choice(correlated_data, size=nsubset, replace=False)

# Report the correlation coefficient of the sample
sample_corr = np.corrcoef(subset[:, 0], subset[:, 1])[0, 1]
print(f'The correlation coefficient of our sample is {sample_corr:.3f}.')

# Plot both the population and our sample
plt.scatter(correlated_data[:, 0], correlated_data[:, 1], marker='x', c='black')
plt.scatter(subset[:, 0], subset[:, 1], c='red')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(['True population', 'Sample'])
plt.show()

The correlation coefficient of our sample is -0.619.

Maintenant, le bootstrap. Chaque rééchantillon tire len(subset) paires de l’échantillon lui-même, avec remise. La taille du rééchantillon égale celle de l’échantillon d’origine : c’est ce qui fait que la dispersion de la statistique rééchantillonnée imite la variabilité d’échantillonnage de l’estimation d’origine.

number_runs = 1000

# Array to record the correlation coefficient of each resample
corr_coef_collector = np.zeros(number_runs)

# The bootstrap resample size equals the original sample size.
length_sub = len(subset)

for i in range(number_runs):
    # Draw length_sub pairs from the sample, WITH REPLACEMENT
    new_pairs = rng.choice(subset, size=length_sub, replace=True)
    corr_coef_collector[i] = np.corrcoef(new_pairs[:, 0], new_pairs[:, 1])[0, 1]

# Plot the bootstrap distribution
plt.hist(corr_coef_collector, color='black')
plt.xlabel('Correlation coefficient')
plt.ylabel('Count')

plt.axvline(x=correlation_matrix[0, 1], color='red')
plt.axvline(x=sample_corr, color='orange')
plt.axvline(x=np.median(corr_coef_collector), color='blue')
plt.legend(['True correlation coefficient', 'Sample correlation coefficient',
            'Median of bootstrap estimates'])
plt.show()

La médiane des estimations bootstrap se situe près du coefficient de corrélation de l’échantillon, pas de la vraie valeur de la population. Le bootstrap ne peut pas corriger le biais d’un petit échantillon ; il quantifie l’incertitude autour de l’estimation dont vous disposez.

Que se passe-t-il si vous augmentez la taille de l’échantillon d’origine ?

Et si la taille du rééchantillon diffère de celle de l’échantillon ?¶

Le bootstrap se définit avec une taille de rééchantillon égale à la taille de l’échantillon. À titre d’expérience, nous enfreignons délibérément cette règle et faisons varier la taille du rééchantillon. Observez la dispersion de la distribution bootstrap.

resample_sizes = [10, 25, 50, 200]  # 50 is the actual sample size

fig, ax = plt.subplots(figsize=(8, 4))
for m in resample_sizes:
    stats_m = np.zeros(number_runs)
    for i in range(number_runs):
        new_pairs = rng.choice(subset, size=m, replace=True)
        stats_m[i] = np.corrcoef(new_pairs[:, 0], new_pairs[:, 1])[0, 1]
    ax.hist(stats_m, bins=30, histtype='step', lw=2,
            label=f'resample size {m}, std {np.std(stats_m):.3f}')

ax.axvline(sample_corr, color='k', ls='--', label='sample correlation')
ax.set_xlabel('Correlation coefficient')
ax.set_ylabel('Count')
ax.legend()
plt.show()

Des rééchantillons plus petits produisent une dispersion plus large de la statistique : une corrélation estimée sur 10 paires est plus bruitée qu’une corrélation estimée sur 50. Des rééchantillons plus grands que l’échantillon produisent une dispersion trop étroite, qui sous-estime l’incertitude réelle. Seule une taille de rééchantillon égale à la taille de l’échantillon d’origine reproduit la variabilité d’échantillonnage de l’estimation que vous avez réellement faite.

1.3 Monte-Carlo¶

Nommées d’après le casino de Monaco, les méthodes de Monte-Carlo consistent à simuler de nouvelles données à partir d’un modèle statistique connu (ou supposé !). Contrairement aux deux exemples précédents, nous ne tirons pas d’un échantillon existant.

Les techniques de Monte-Carlo ont de nombreuses applications : évaluation probabiliste des risques, propagation d’incertitude dans les modèles, évaluation de systèmes trop compliqués pour une analyse en forme close.

Ici, nous illustrons l’échantillonnage de Monte-Carlo en estimant π.

Nous partons d’une idée centrale : le rapport de l’aire d’un cercle à l’aire du carré qui le borne vaut $\pi/4$ .

Nous imaginons alors un cercle de rayon 1 inscrit dans un carré dont les côtés vont de -1 à 1.

# We draw (x, y) points from a *uniform* distribution and determine whether each point
# sits within the circle. Every point lands in the square; only some land in the circle.
in_circle = np.empty([0, 2])
in_square = np.empty([0, 2])

# Generate the samples. Keep the number small at first.
number_runs = 50

for _ in range(number_runs):  # note that _ avoids creating a loop variable we never use
    x = rng.uniform(low=-1, high=1)
    y = rng.uniform(low=-1, high=1)
    # How far is this point from the origin?
    origin_dist = x**2 + y**2
    # If origin_dist is less than 1, the point is inside the circle
    if origin_dist <= 1:
        in_circle = np.append(in_circle, [[x, y]], axis=0)
    in_square = np.append(in_square, [[x, y]], axis=0)

# Visualize what we just did
plt.scatter(in_square[:, 0], in_square[:, 1], marker='x', c='black')
plt.scatter(in_circle[:, 0], in_circle[:, 1], marker='o', c='red')
plt.xlabel('X coordinate')
plt.ylabel('Y coordinate')
plt.legend(['Points in square', 'Points in circle'], loc='center left', bbox_to_anchor=(1, 0.5))
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()

Nous estimons maintenant π en utilisant les comptes de points dans le cercle et dans le carré comme approximations de leurs aires.

pi_est = 4 * (len(in_circle) / len(in_square))
print(f'We estimate the value of pi to be: {pi_est}.')

We estimate the value of pi to be: 3.2.

Avec quelques tirages, nous n’obtenons pas une bonne réponse.

Une approche Monte-Carlo demande plus de tirages pour converger.

Explorez comment le nombre de tirages change votre estimation de π (et comment votre calcul converge).

2. Le rééchantillonnage pour une inférence de modèle robuste (niveau 2)¶

Le plan : ajuster une tendance linéaire à des données de position GNSS et utiliser le bootstrap pour attribuer une incertitude à la pente (la vitesse de la plaque). Nous commençons par une série synthétique, dont la vraie vitesse est connue, afin de comparer la distribution bootstrap à la bonne réponse. Puis nous répétons l’analyse sur des données réelles de la station P395, dans le Nord-Ouest Pacifique.

2.1 Série GNSS synthétique à vitesse connue¶

L’utilitaire mlgeo_synth.gnss_series construit une série synthétique de déplacements GNSS journaliers avec une vérité terrain connue : une tendance linéaire, des cycles saisonniers et un bruit réaliste (bruit blanc, bruit de scintillation (flicker) et marche aléatoire). Nous fixons la vraie vitesse à 12 mm/an.

true_velocity = 12.0  # mm/yr, our ground truth

gnss = mlgeo_synth.gnss_series(n_years=10, velocity_mm_yr=true_velocity, seed=42)
gnss.head()

# Time in years since the first sample
t_syn = (gnss['date'] - gnss['date'].iloc[0]).dt.days / 365.25
d_syn = gnss['disp_mm']

plt.plot(t_syn, d_syn, lw=0.5, label='synthetic displacement')
plt.plot(t_syn, gnss['trend_mm'], 'r', label='true trend (12 mm/yr)')
plt.xlabel('Time (years)')
plt.ylabel('Displacement (mm)')
plt.legend()
plt.show()

Ajustez une droite avec scipy.stats.linregress. La pente est notre estimation de la vitesse.

from scipy import stats

fit = stats.linregress(t_syn, d_syn)
print(f'Estimated velocity: {fit.slope:.3f} mm/yr (true value: {true_velocity} mm/yr)')

Estimated velocity: 11.838 mm/yr (true value: 12.0 mm/yr)

L’estimation ponctuelle est proche de la vérité, mais quelle confiance lui accorder ? Bootstrap : rééchantillonnez les paires (temps, déplacement) avec remise — taille de rééchantillon égale à la longueur des données —, réajustez la droite à chaque fois et collectez les pentes.

k = 1000
n_syn = len(t_syn)
t_arr = t_syn.to_numpy()
d_arr = d_syn.to_numpy()

vel_syn = np.zeros(k)
for j in range(k):
    ii = rng.integers(0, n_syn, size=n_syn)  # indices drawn with replacement
    vel_syn[j] = stats.linregress(t_arr[ii], d_arr[ii]).slope

print(f'Bootstrap mean velocity: {np.mean(vel_syn):.3f} mm/yr, '
      f'standard deviation: {np.std(vel_syn):.3f} mm/yr')

plt.hist(vel_syn, bins=30, color='black')
plt.axvline(true_velocity, color='red', label='true velocity')
plt.axvline(np.mean(vel_syn), color='orange', label='bootstrap mean')
plt.xlabel('Velocity (mm/yr)')
plt.ylabel('Count')
plt.legend()
plt.show()

Bootstrap mean velocity: 11.838 mm/yr, standard deviation: 0.019 mm/yr

Deux choses à remarquer.

D’abord, la distribution bootstrap est centrée sur la pente estimée, pas sur la vérité. Le bootstrap quantifie la variabilité de l’estimateur ; il ne peut pas supprimer l’écart entre l’estimation et la vraie vitesse que cette réalisation particulière du bruit a produit.

Ensuite, la distribution est très étroite — et la vraie vitesse se trouve à plusieurs écarts-types bootstrap de son centre. La barre d’erreur est trop petite. La raison : rééchantillonner les paires traite les résidus comme indépendants, alors que le bruit GNSS est corrélé dans le temps (scintillation et marche aléatoire). Le bootstrap par paires sous-estime donc l’incertitude réelle. Des schémas plus avancés (bootstrap par blocs) rééchantillonnent des tronçons contigus de la série pour préserver la corrélation. C’est la vérité terrain qui a révélé le problème : avec les seules données réelles, l’histogramme étroit aurait paru rassurant.

2.2 Mouvement des plaques à partir de données géodésiques réelles¶

Place au réel. Nous utilisons une série temporelle GNSS de la station P395, dans le Nord-Ouest Pacifique, et estimons le mouvement à long terme dû à la zone de subduction de Cascadia.

Nous téléchargeons la série depuis le centre de données de l’Université du Nevada à Reno. Le fichier tenv3 est délimité par des espaces, avec une ligne d’en-tête. Les colonnes utiles sont l’année décimale (yyyy.yyyy) et les positions est, nord et verticale en mètres (__east(m), _north(m), ____up(m)).

sta = "P395"
url = f"https://geodesy.unr.edu/gps_timeseries/IGS20/tenv3/IGS20/{sta}.tenv3"
print(url)

os.makedirs('data', exist_ok=True)
fname = f'data/{sta}.tenv3'
r = requests.get(url, timeout=60)
r.raise_for_status()
with open(fname, 'wb') as f:
    f.write(r.content)

# Whitespace-delimited file; the first line is the header.
df = pd.read_csv(fname, sep=r'\s+')
df.head()

https://geodesy.unr.edu/gps_timeseries/IGS20/tenv3/IGS20/P395.tenv3

# Keep only the columns we need and give them simpler names.
df = df[['yyyy.yyyy', '__east(m)', '_north(m)', '____up(m)']].rename(
    columns={'yyyy.yyyy': 'decimal year',
             '__east(m)': 'delta e (m)',
             '_north(m)': 'delta n (m)',
             '____up(m)': 'delta v (m)'})

# Drop rows with missing values. dropna returns a new frame: assign the result.
df = df.dropna()
df.head()

# Reference each component to the first epoch so positions start at zero.
df['new delta e (m)'] = df['delta e (m)'] - df['delta e (m)'].values[0]
df['new delta n (m)'] = df['delta n (m)'] - df['delta n (m)'].values[0]
df['new delta v (m)'] = df['delta v (m)'] - df['delta v (m)'].values[0]
df.head()

plt.plot(df['decimal year'], df['new delta e (m)'], label='East displacement')
plt.xlabel('Year')
plt.ylabel('Displacement (m)')
plt.legend()
plt.show()

2.3 Régression linéaire¶

Il y a une tendance linéaire nette dans les positions horizontales. Nous pouvons ajuster les données avec :

$E(t) = V_e t + u_e,$
$N(t) = V_n t + u_n,$

où $t$ est le temps. Nous régressons les données pour trouver les coefficients $V_e$ , $u_e$ , $V_n$ , $u_n$ . Les déplacements sont surtout vers l’ouest, donc nous nous concentrons sur la composante Est $E(t)$ pour cet exercice. Les coefficients $u_e$ et $u_n$ sont les ordonnées à l’origine en $t=0$ . Ils ne sont pas nuls ici parce que $t$ commence en 2006. Les coefficients $V_e$ et $V_n$ ont la dimension de vitesses :

$V_e \sim E(t) / t$ , $V_n \sim N(t) / t$ ,

cet exemple nous permet donc de discuter une régression linéaire simple et le rééchantillonnage. Nous utilisons à la fois une fonction SciPy et une fonction scikit-learn.

Pour mesurer la performance de l’ajustement, nous mesurons dans quelle mesure la variance est réduite en ajustant les données (les points du nuage) au modèle. La variance est :

$\text{Var}(x) = 1/n \sum_{i=1}^n (x_i-\hat{x})^2$ ,

où $\hat{x}$ est la moyenne de $x$ . En ajustant la régression, nous prédisons les valeurs $x_{pred}$ . Les résidus sont les différences entre les données et les valeurs prédites : $e = x - x_{pred}$ . Le $R^2$ ou coefficient de détermination est :

$R^2 = 1 - \text{Var}(x-x_{pred}) /\text{Var}(x) = 1 - \text{Var}(e) /\text{Var}(x)$

Plus l’erreur est petite, « meilleur » est l’ajustement (nous discuterons plus loin du fait qu’un ajustement peut être trop bon !), et plus $R^2$ est proche de un.

# Linear regression: displacement = velocity * time + intercept, East component.
Ve, intercept, r_value, p_value, std_err = stats.linregress(df['decimal year'],
                                                            df['new delta e (m)'])
print(sta, "overall plate motion there", Ve, 'm/year')
print("parameters: correlation coefficient %4.2f, P-value %4.2f, standard error of the slope %g"
      % (r_value, p_value, std_err))

P395 overall plate motion there -0.006540884541421708 m/year
parameters: correlation coefficient -1.00, P-value 0.00, standard error of the slope 5.78311e-06

Nous pouvons aussi utiliser le paquet scikit-learn :

from sklearn.linear_model import LinearRegression

# Convert the data into numpy arrays. Reshaping to (n, 1) is required by scikit-learn.
E = np.asarray(df['new delta e (m)']).reshape(-1, 1)
t = np.asarray(df['decimal year']).reshape(-1, 1)

# Perform the linear regression on the entire available data
regr = LinearRegression()
regr.fit(t, E)
Epred = regr.predict(t)

# The coefficients
print('Coefficient / velocity eastward (m/year): ', regr.coef_[0][0])

# Plot the data and the fit
plt.plot(t, E, 'b', label='data')
plt.plot(t, Epred, 'r', label='linear fit')
plt.xlabel('Year')
plt.ylabel('East displacement (m)')
plt.legend()
plt.show()

Coefficient / velocity eastward (m/year):  -0.00654088454142171

Pour évaluer les erreurs de l’ajustement du modèle avec sklearn, nous utilisons les fonctions suivantes :

from sklearn.metrics import mean_squared_error, r2_score

# The mean squared error
print('Mean squared error (m^2): %.6f' % mean_squared_error(E, Epred))
# The coefficient of determination: 1 is the perfect prediction
print('Coefficient of determination: %.2f' % r2_score(E, Epred))

Mean squared error (m^2): 0.000009
Coefficient of determination: 0.99

2.4 Bootstrap de la vitesse¶

Nous utilisons maintenant le bootstrap pour estimer la pente de la régression sur de nombreux jeux de données rééchantillonnés, exactement comme pour la série synthétique.

Scikit-learn fournit resample dans le module utils. Veillez à utiliser replace=True et une taille de rééchantillon égale à la longueur des données (le défaut). Pour des résultats reproductibles, vous pouvez passer un random_state fixe. Le bootstrap se répète en général de nombreuses fois (au contraire de la validation croisée en K plis (K-fold), le schéma d’évaluation de modèles du chapitre 3.8, qui scinde les données en un nombre fixe de plis sans recouvrement).

from sklearn.utils import resample

k = 1000
vel = np.zeros(k)  # initialize a vector to store the regression slopes
for i in range(k):
    ii = resample(np.arange(len(E)), replace=True, n_samples=len(E),
                  random_state=i)  # new indices
    E_b, t_b = E[ii], t[ii]
    # Fit the resampled data
    regr = LinearRegression()
    regr.fit(t_b, E_b)
    vel[i] = regr.coef_[0][0]

# The data shows a clear trend, so the slope estimates are close to each other:
print("mean of the velocity estimates %g m/yr and standard deviation %g m/yr"
      % (np.mean(vel), np.std(vel)))

plt.hist(vel, 10)
plt.title('Distribution of eastward velocities (m/year)')
plt.xlabel('Velocity (m/year)')
plt.ylabel('Count')
plt.grid(True)
plt.show()

mean of the velocity estimates -0.00654102 m/yr and standard deviation 5.34869e-06 m/yr

La dispersion bootstrap sur les données réelles est petite parce que la tendance domine le bruit, exactement comme dans le cas synthétique. La même réserve s’applique : le bruit GNSS est corrélé dans le temps, donc traitez cette barre d’erreur comme une borne inférieure.

3. Rééchantillonnage de signaux et données irrégulières (niveau 3)¶

Les sections 1 et 2 utilisaient « rééchantillonnage » au sens statistique : retirer d’un échantillon pour mesurer une incertitude. Le mot a un second sens que tout flux de capteurs finit tôt ou tard par vous imposer : changer l’échantillonnage d’une série temporelle — sous-échantillonner un enregistrement à haute cadence, combler ou refuser de combler des lacunes, et placer des observations irrégulières sur une grille régulière. Les deux sens partagent un même piège : appliqués naïvement, ils fabriquent un signal qui n’a jamais été mesuré.

Cette section travaille sur trois flux de données qui couvrent l’essentiel de ce que vous rencontrerez en pratique :

une série régulière à haute cadence (un marégraphe horaire) que nous sous-échantillonnons en valeurs journalières ;
une série régulière avec interruptions (un enregistrement GNSS journalier à lacunes) que nous interpolons selon une politique de lacunes explicite ;
un flux de points épars et irrégulier (un réseau de puits d’eaux souterraines suivi sur plusieurs décennies) où rien dans l’échantillonnage n’est régulier et où les choix d’agrégation dominent le résultat.

Chacun est synthétique avec une vérité terrain connue, donc chaque réparation est notée. Nous terminons en revenant sur la barre d’erreur bootstrap trop étroite de la section 2.1, que nous corrigeons avec un bootstrap par blocs.

3.1 Sous-échantillonnage et repliement¶

Le sous-échantillonnage a l’air inoffensif : garder un échantillon sur $q$ , jeter le reste. Il ne l’est pas. Une série échantillonnée à l’intervalle $\Delta t$ ne peut représenter que des fréquences jusqu’à la fréquence de Nyquist $f_N = 1/(2\Delta t)$ . Tout signal au-dessus de la nouvelle fréquence de Nyquist ne disparaît pas quand vous sous-échantillonnez — il se replie vers une fréquence plus basse, en se faisant passer pour un signal qui n’a jamais existé.

Le banc d’essai idéal est un marégraphe — celui de Brest, l’une des plus longues séries marégraphiques au monde, en est l’archétype. mlgeo_synth.tide_gauge_series génère un niveau de la mer horaire à partir de quatre composantes astronomiques, plus une tendance, un cycle saisonnier et un bruit météorologique — et renvoie les composantes de marée comme vérité terrain. Supposons que nous voulions une série journalière du niveau de la mer pour étudier le signal lent (subtidal).

# Six months of hourly sea level
tide, tide_truth = mlgeo_synth.tide_gauge_series(n_days=180, seed=42)
sl = tide.set_index('time')['sea_level_m']

fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(sl.iloc[:24 * 10], lw=0.8)
ax.set_ylabel('sea level (m)')
ax.set_title('Hourly tide-gauge record, first 10 days')
plt.tight_layout()
plt.show()

tide_truth['constituents']

La composante dominante est M2, la marée semi-diurne lunaire principale, de période 12,42 heures — une fréquence de 1,93 cycle par jour. Une série journalière a une fréquence de Nyquist de 0,5 cycle par jour, donc la marée entière vit au-dessus de la nouvelle fréquence de Nyquist. Si nous gardons un échantillon par jour (la lecture de minuit, par exemple), M2 se replie à $|1.93 - 2| = 0.068$ cycle par jour : une oscillation parasite de période 14,8 jours, avec la pleine amplitude de marée d’environ 0,8 m.

La correction est la règle que tout sous-échantillonnage doit suivre : d’abord filtrer passe-bas sous la nouvelle fréquence de Nyquist, puis sous-échantillonner. Une moyenne journalière est un filtre passe-bas grossier (une fenêtre rectangulaire (boxcar) de 24 heures) et supprime déjà l’essentiel de la marée ; scipy.signal.decimate applique un vrai filtre anti-repliement avant de sous-échantillonner. Nous notons les trois versions contre le vrai niveau de la mer journalier sans marée, que nous pouvons calculer exactement puisque le générateur a renvoyé la marée dans une colonne séparée.

from scipy import signal

# Ground truth: the daily mean of the tide-free sea level
subtidal = (tide.set_index('time')['sea_level_m']
            - tide.set_index('time')['tide_m']).resample('D').mean()

# WRONG: keep one sample per day (midnight), discard the rest
naive = sl.iloc[::24]

# Crude anti-alias: daily mean (a 24-h boxcar low-pass, then subsample)
daily_mean = sl.resample('D').mean()

# Proper anti-alias: decimate in two stages (scipy recommends factors <= 13)
dec = signal.decimate(signal.decimate(sl.to_numpy(), 4, ftype='fir', zero_phase=True),
                      6, ftype='fir', zero_phase=True)
dec = pd.Series(dec, index=sl.index[::24])

fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(naive, color='tab:red', lw=1, label='midnight sample (aliased)')
ax.plot(daily_mean, color='tab:orange', lw=1.2, label='daily mean')
ax.plot(dec, color='tab:blue', lw=1.2, label='decimate (anti-aliased)')
ax.plot(subtidal, 'k--', lw=1.2, label='true subtidal sea level')
ax.set_ylabel('sea level (m)')
ax.legend(ncols=2)
ax.set_title('Three ways to make a daily series from an hourly one')
plt.tight_layout()
plt.show()

for name, s_daily in [('midnight sample', naive), ('daily mean', daily_mean),
                      ('decimate', dec)]:
    err = (s_daily - subtidal).dropna()
    print(f'{name:16s} RMS error vs true subtidal signal: {np.sqrt((err**2).mean()):.3f} m')

midnight sample  RMS error vs true subtidal signal: 0.589 m
daily mean       RMS error vs true subtidal signal: 0.020 m
decimate         RMS error vs true subtidal signal: 0.017 m

Les échantillons de minuit portent une oscillation d’environ 15 jours de plus d’un demi-mètre qui n’existe pas dans l’océan subtidal — c’est la marée M2 repliée, et son erreur RMS est trente fois plus grande que celle des deux versions anti-repliées. Rien dans la série naïve n’a l’air faux ; l’ondulation bimensuelle ressemble même à un signal océanique plausible. C’est ce qui rend le repliement dangereux : l’artefact est physiquement habillé. L’altimétrie satellitaire vit exactement avec ce problème — l’orbite des missions franco-américaines TOPEX/Jason échantillonne chaque point tous les ~10 jours, repliant M2 en un signal de 62 jours qu’il faut modéliser et retirer.

pandas.DataFrame.resample('D').mean() — la ligne unique vers laquelle vous tendrez le plus souvent — est déjà un filtre anti-repliement convenable pour cet usage. La règle à intérioriser : avant de réduire la fréquence d’échantillonnage, demandez-vous ce qui vit au-dessus de la nouvelle fréquence de Nyquist et retirez-le. Si la réponse est « rien », dites-le explicitement dans votre fiche de données.

3.2 Lacunes : l’interpolation est une décision, pas un défaut¶

Les vrais flux journaliers arrivent avec des interruptions. mlgeo_synth.degrade_series injecte des lacunes dans la série GNSS synthétique de la section 2.1 — plus une interruption de 150 jours que nous plaçons, délibérément, sur un séisme (une station mise hors service par la secousse qu’elle devait enregistrer n’a rien d’hypothétique). La fonction renvoie la vérité non censurée, pour que nous puissions noter toute réparation.

# The same 10-yr, 12 mm/yr station, now with a coseismic step at day 2000
gnss_eq = mlgeo_synth.gnss_series(n_years=10, velocity_mm_yr=12.0,
                                  eq_day=2000, coseismic_mm=25.0, seed=42)

# Degrade it: an imposed 150-day outage swallowing the earthquake, plus 8 random gaps
broken, truth = mlgeo_synth.degrade_series(gnss_eq, gap_windows=[(1950, 2100)],
                                           n_random_gaps=8, gap_days=(2, 25), seed=13)

s = broken.set_index('date')['disp_mm']
clean = pd.Series(truth['clean'].to_numpy(), index=s.index)

fig, ax = plt.subplots(figsize=(11, 3.5))
ax.plot(s, lw=0.5, label='observed (gaps are blank)')
for s0, e0 in truth['gap_windows']:
    ax.axvspan(s.index[s0], s.index[e0 - 1], color='tab:red', alpha=0.15)
ax.set_ylabel('displacement (mm)')
ax.legend()
ax.set_title(f'Degraded GNSS series: {len(truth["gap_windows"])} gaps (shaded)')
plt.tight_layout()
plt.show()

La ligne unique tentante est s.interpolate() : relier les points à travers chaque lacune. Avant de lui faire confiance, mesurez ce qu’elle coûte — interpolez tout et comparez à la vérité, lacune par lacune.

filled_all = s.interpolate(method='time')

noise_std = (gnss_eq['disp_mm'] - gnss_eq['trend_mm'] - gnss_eq['seasonal_mm']
             - gnss_eq['eq_mm']).std()
print(f'daily noise level of this series: {noise_std:.2f} mm RMS\n')

print('gap length   RMS error of linear interpolation inside the gap')
for s0, e0 in sorted(truth['gap_windows'], key=lambda w: w[1] - w[0]):
    idx = s.index[s0:e0]
    rms = np.sqrt(((filled_all - clean)[idx] ** 2).mean())
    print(f'{e0 - s0:7d} d   {rms:5.2f} mm')

daily noise level of this series: 2.41 mm RMS

gap length   RMS error of linear interpolation inside the gap
      3 d    0.74 mm
      6 d    1.84 mm
     16 d    1.60 mm
     20 d    1.35 mm
     21 d    1.87 mm
     21 d    1.79 mm
     23 d    1.70 mm
     24 d    2.37 mm
    150 d    8.09 mm

Les lacunes courtes s’interpolent au niveau du bruit ou en dessous : en quelques jours, la tendance et le cycle saisonnier bougent à peine, donc une droite vaut les données. La lacune de 150 jours est différente — son erreur dépasse trois fois le bruit, parce que l’interpolation a tracé une rampe lisse à travers un saut cosismique de 25 mm qu’elle n’avait aucun moyen de connaître. Les valeurs comblées ne sont pas bruitées ; elles sont fausses avec assurance, et tout code en aval (un ajustement de tendance, une fenêtre de caractéristiques pour le ML) les traitera comme des mesures.

Nous énonçons donc une politique de lacunes comme une décision explicite plutôt que comme un défaut de bibliothèque :

Interpoler les lacunes de 10 jours ou moins. Justification, d’après le tableau ci-dessus : à 12 mm/an et avec une amplitude saisonnière d’environ 3 mm, le signal déterministe bouge bien moins que le niveau de bruit de 2,4 mm en 10 jours, donc l’erreur d’interpolation est bornée par le bruit. Le seuil est fixé par les taux de signal de cette station — une station au mouvement plus rapide ou aux oscillations saisonnières plus fortes mérite un seuil plus court, et le seuil a sa place dans la fiche de données.
Masquer les lacunes plus longues comme manquantes. Un NaN est une déclaration honnête : nous ne savons pas ce qui s’est passé là — et dans cet enregistrement, il s’est bel et bien passé quelque chose.

max_gap_days = 10  # the decision, justified above

# Length of the gap each missing sample belongs to
isna = s.isna()
gap_id = (isna != isna.shift()).cumsum()
gap_len = isna.groupby(gap_id).transform('sum').where(isna, 0)

# Fill short gaps, keep long gaps as NaN
repaired = filled_all.where(~(isna & (gap_len > max_gap_days)))

short_filled = isna & (gap_len <= max_gap_days)
long_masked = isna & (gap_len > max_gap_days)
rms_short = np.sqrt(((repaired - clean)[short_filled] ** 2).mean())
rms_if_filled = np.sqrt(((filled_all - clean)[long_masked] ** 2).mean())

print(f'samples filled (short gaps): {short_filled.sum()},  RMS error {rms_short:.2f} mm '
      f'(noise level {noise_std:.2f} mm)')
print(f'samples masked (long gaps):  {long_masked.sum()},  RMS error if we had '
      f'interpolated them: {rms_if_filled:.2f} mm')

# Zoom on the long gap: what interpolation would have fabricated
zoom = slice('2020-01-01', '2020-12-31')
fig, ax = plt.subplots(figsize=(10, 3.5))
ax.plot(clean[zoom], color='gray', lw=0.6, label='truth (never observed)')
ax.plot(filled_all[zoom].where(long_masked[zoom]), 'r--', lw=1.5,
        label='linear interpolation (fabricated)')
ax.plot(repaired[zoom], lw=0.8, color='tab:blue', label='policy: filled + masked')
ax.set_ylabel('displacement (mm)')
ax.legend()
ax.set_title('The long gap hides an earthquake; interpolation invents a smooth story')
plt.tight_layout()
plt.show()

samples filled (short gaps): 9,  RMS error 1.56 mm (noise level 2.41 mm)
samples masked (long gaps):  275,  RMS error if we had interpolated them: 6.10 mm

Le résultat noté : l’interpolation des lacunes courtes coûte moins que le bruit, et le masque refuse d’inventer les 150 jours que nous n’avons jamais vus. Quand cette série alimentera plus tard un modèle, les NaN forceront un choix documenté (écarter la fenêtre, la signaler, imputer avec une incertitude) au lieu de propager silencieusement de la fiction. C’est le motif de chaque réparation de cette section : corriger ce que les données contraignent, masquer ce qu’elles ne contraignent pas, et écrire le seuil noir sur blanc.

3.3 Points épars irréguliers : le réseau multi-puits¶

Le troisième flux n’a pas de grille de départ. mlgeo_synth.well_table imite quarante ans de mesures de niveau d’eau sur 25 puits de suivi : tous les puits partagent un même signal régional (une recharge saisonnière superposée à un lent déclin), mais chacun a son propre décalage de référence (des mètres !), son propre niveau de bruit, sa propre période d’activité, une lacune pluriannuelle et des dates de visite inégales. Quelques puits n’ont été visités qu’une poignée de fois. Voilà à quoi ressemblent réellement l’hydrologie opérationnelle, le suivi géotechnique et les archives héritées — le réseau piézométrique national français, dont les données sont diffusées par la banque ADES, en est un exemple grandeur nature.

wells, wtruth = mlgeo_synth.well_table(seed=42)
print(f"{wells['well_id'].nunique()} wells, {len(wells)} observations, "
      f"{wells['date'].min():%Y} to {wells['date'].max():%Y}")

fig, ax = plt.subplots(figsize=(11, 4))
sc = ax.scatter(wells['date'], wells['head_m'], c=wells['well_id'], s=4, cmap='tab20')
ax.set_ylabel('head (m, arbitrary regional datum)')
ax.set_title('25 wells: shared regional signal buried under per-well offsets and gaps')
plt.tight_layout()
plt.show()

25 wells, 2902 observations, 1981 to 2019

L’objectif : retrouver le signal de charge hydraulique régional — la tendance et le cycle saisonnier partagés — en série trimestrielle. Le geste naïf est resample('QS').mean() sur toutes les observations. Regardez ce qu’il fait : chaque fois qu’un puits à fort décalage de référence entre dans l’enregistrement ou en sort (et ils entrent et sortent tous à des moments différents), la moyenne saute d’une fraction de ce décalage. Le « signal régional » qu’il produit est surtout une histoire de quels puits étaient visités.

La version consciente des lacunes fait deux gestes :

Travailler en anomalies. Soustraire d’abord la moyenne propre de chaque puits, pour que les décalages de référence métriques s’annulent avant toute moyenne. (C’est exactement ainsi que les séries de température globale sont construites à partir des stations météo.)
Pondérer par la qualité de mesure. Chaque observation porte un sigma_m déclaré ; pondérer par $1/\sigma^2$ — le poids par inverse de la variance, qui minimise la variance de l’estimation combinée — empêche une lecture au ruban gradué de diluer un enregistrement de capteur de pression.

Les anomalies laissent encore un petit biais — chaque puits échantillonne un morceau différent du déclin sur 40 ans, donc sa propre moyenne absorbe un segment de tendance légèrement différent — mais ce résidu se compte en décimètres, pas dans les mètres que les décalages injecteraient. Nous notons les deux contre truth["regional"], à une constante près (la référence est arbitraire, donc nous comparons toutes les séries après retrait de leur moyenne).

# Naive: average all raw heads in each quarter
naive_q = wells.set_index('date')['head_m'].resample('QS').mean()

# Gap-aware: per-well anomalies, inverse-variance weights, quarterly aggregation
w = wells.assign(
    anom_m=wells['head_m'] - wells.groupby('well_id')['head_m'].transform('mean'),
    weight=1.0 / wells['sigma_m'] ** 2,
    quarter=wells['date'].dt.to_period('Q').dt.start_time,
)
w['wx'] = w['weight'] * w['anom_m']
g = w.groupby('quarter')
aware_q = g['wx'].sum() / g['weight'].sum()

# Grade against the true regional signal (all series demeaned: the datum is arbitrary)
def rms_vs_regional(series):
    t_yr = (series.index - wells['date'].min()).days / 365.25
    reg = wtruth['regional'](t_yr)
    return np.sqrt(np.mean(((series - series.mean()) - (reg - reg.mean())) ** 2))

t_yr = (aware_q.index - wells['date'].min()).days / 365.25
regional_true = pd.Series(wtruth['regional'](t_yr), index=aware_q.index)

fig, ax = plt.subplots(2, 1, figsize=(11, 6), sharex=True)
ax[0].plot(naive_q - naive_q.mean(), color='tab:red', lw=0.8)
ax[0].set_title(f'Naive quarterly mean — RMS error {rms_vs_regional(naive_q):.2f} m')
ax[1].plot(aware_q - aware_q.mean(), color='tab:blue', lw=0.8,
           label=r'anomaly + 1/$\sigma^2$ aggregation')
ax[1].plot(regional_true - regional_true.mean(), 'k--', lw=1,
           label='true regional signal')
ax[1].set_title(f'Gap-aware aggregation — RMS error {rms_vs_regional(aware_q):.2f} m')
ax[1].legend()
for a in ax:
    a.set_ylabel('head anomaly (m)')
plt.tight_layout()
plt.show()

La moyenne naïve se trompe d’un facteur de plusieurs unités — des sauts métriques produits entièrement par les entrées et sorties de puits dans l’enregistrement — tandis que la série pondérée en anomalies suit le vrai déclin régional et son cycle saisonnier à environ un demi-mètre RMS près — surtout le biais résiduel de segments de tendance noté ci-dessus. Rien de sophistiqué ne s’est produit : toute l’amélioration vient du refus de moyenner des choses incompatibles. Avant d’agréger tout flux multi-sites irrégulier, demandez-vous ce qui entre dans la moyenne et en sort quand la composition change, et retirez d’abord les niveaux propres à chaque site.

3.4 Boucler la boucle : le bootstrap par blocs mobiles¶

La section 2.1 s’est terminée sur un avertissement : le bootstrap par paires donnait une barre d’erreur de vitesse si étroite que la vraie vitesse se trouvait loin en dehors, parce que rééchantillonner des jours individuels détruit la corrélation temporelle du bruit GNSS, et le bruit corrélé est exactement ce qui rend une tendance incertaine. La correction est le bootstrap par blocs mobiles : au lieu de rééchantillonner des jours, on rééchantillonne des blocs contigus de résidus, de sorte que chaque rééchantillon préserve la corrélation du bruit jusqu’à la longueur du bloc. Nous ajustons la droite une fois, rééchantillonnons des blocs de ses résidus, les rajoutons à la droite ajustée et réajustons.

La longueur de bloc est — encore une fois — une décision énoncée : elle doit dépasser le temps de corrélation du bruit qui vous importe. Nous utilisons 100 jours, confortablement plus long que la corrélation du bruit de scintillation aux périodes qui comptent pour une tendance décennale ; vous pouvez vérifier la sensibilité en relançant avec 50 ou 200.

block = 100                     # days per block: the decision
n_blocks = int(np.ceil(n_syn / block))
line = fit.intercept + fit.slope * t_arr
resid = d_arr - line

vel_block = np.zeros(k)
for j in range(k):
    starts = rng.integers(0, n_syn - block, size=n_blocks)
    boot_resid = np.concatenate([resid[s0:s0 + block] for s0 in starts])[:n_syn]
    vel_block[j] = stats.linregress(t_arr, line + boot_resid).slope

print(f'pair bootstrap:  std {np.std(vel_syn):.3f} mm/yr, '
      f'true velocity sits {abs(fit.slope - true_velocity) / np.std(vel_syn):.1f} sigma out')
print(f'block bootstrap: std {np.std(vel_block):.3f} mm/yr, '
      f'true velocity sits {abs(fit.slope - true_velocity) / np.std(vel_block):.1f} sigma out')

fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(vel_syn, bins=30, color='black', alpha=0.7, label='pair bootstrap (sec. 2.1)')
ax.hist(vel_block, bins=30, color='tab:blue', alpha=0.6, label=f'block bootstrap ({block}-day blocks)')
ax.axvline(true_velocity, color='red', label='true velocity')
ax.set_xlabel('velocity (mm/yr)')
ax.set_ylabel('count')
ax.legend()
plt.show()

pair bootstrap:  std 0.019 mm/yr, true velocity sits 8.6 sigma out
block bootstrap: std 0.148 mm/yr, true velocity sits 1.1 sigma out

Le bootstrap par blocs élargit la barre d’erreur d’un facteur d’environ sept — et la vraie vitesse se trouve maintenant à environ un écart-type de l’estimation, ce à quoi ressemble une barre d’erreur honnête. Rien n’a changé dans les données ; seul le rééchantillonnage a respecté la corrélation que le bruit possède réellement. L’histogramme étroit du bootstrap par paires n’était pas prudent, il était faux — et sur des données réelles, sans vérité terrain pour le signaler, il aurait été publié.

La même réserve se transfère maintenant à l’enregistrement réel P395 de la section 2.4 : sa barre d’erreur de bootstrap par paires est une borne inférieure, et un bootstrap par blocs sur ses résidus est l’exercice de prolongement.

4. Exercice¶

Chasse au repliement. Régénérez le marégraphe avec n_days=365 et sous-échantillonnez en gardant l’échantillon de midi au lieu de celui de minuit. La période repliée change-t-elle ? Expliquez pourquoi (ou pourquoi pas) à partir de la formule de repliement.
Sensibilité de la politique de lacunes. Relancez la section 3.2 avec max_gap_days de 3, puis de 30. Rapportez l’erreur RMS et le nombre d’échantillons comblés pour chaque cas. Où placeriez-vous le seuil pour une station se déplaçant à 50 mm/an, et pourquoi ?
Pondération des puits. Dans la section 3.3, supprimez les poids $1/\sigma^2$ (moyenne simple des anomalies). Quelle part de l’amélioration par rapport à la moyenne naïve survit ? Que cela vous dit-il sur celui des deux gestes (anomalies, poids) qui porte la charge pour ce réseau ?
Longueur de bloc. Relancez le bootstrap par blocs avec des longueurs de bloc de 10, 50, 200 et 500 jours et tracez l’écart-type bootstrap en fonction de la longueur de bloc. Expliquez la tendance aux deux extrêmes.