2.7 Considérations statistiques sur les données géoscientifiques et le bruit

Distributions canoniques des données géoscientifiques¶

Distribution normale : utilisée pour des variables comme la température, les variations du niveau de la mer ou la vitesse du vent, dont les valeurs se distribuent symétriquement autour d’une moyenne.
Distribution log-normale : observée dans les phénomènes dont les valeurs ne peuvent pas être négatives et présentent de longues queues, comme l’intensité des pluies, le débit des rivières, la taille des grains et la perméabilité.
Distribution exponentielle : employée pour modéliser les intervalles de temps entre événements, comme le temps entre séismes. Les magnitudes des séismes suivent une forme exponentielle apparentée, la loi de Gutenberg-Richter, que nous détaillons plus bas.
Distribution en loi puissance : présente dans les événements rares de grande ampleur comme les glissements de terrain et les feux de forêt, où les petits événements sont fréquents et les grands événements rares.

🖥️ Diapositives du cours — Séance 07 (mer. 14 oct.)

1. Caractéristiques statistiques¶

Soit $P(z)$ la distribution des données $z$ .

La moyenne¶

Image tirée de ce blog.

La moyenne est la somme des valeurs divisée par le nombre de points de données. C’est le premier moment brut d’une distribution. $\mu = \int_{-\infty}^\infty zP(z)dz$ , où $z$ est la valeur des données (la classe) et $P(z)$ la distribution des données.

La variance¶

La variance est le deuxième moment centré. Centré signifie que la distribution est ramenée autour de la moyenne. Elle calcule l’étalement d’une distribution.

$\sigma^2 = \int_{-\infty}^\infty (z-\mu)^2P(z)dz$

L’écart-type est la racine carrée de la variance, σ. Une variance élevée indique une distribution large.

L’asymétrie (skewness)¶

L’asymétrie est le troisième moment standardisé. Le moment standardisé est mis à l’échelle par l’écart-type. Elle mesure la taille relative des deux queues de la distribution.

$m_3= \int_{-\infty}^\infty \frac{(z - \mu)^3}{\sigma^3}P(z)dz$

Avec l’exposant cubique, l’asymétrie peut être négative.

Image tirée de ce blog.

Une distribution à asymétrie positive est une distribution dont l’essentiel du poids se trouve à la fin de la distribution. Une distribution à asymétrie négative est une distribution dont l’essentiel du poids se trouve au début de la distribution.

Le kurtosis¶

Le kurtosis mesure la taille combinée des deux queues par rapport à l’ensemble de la distribution. C’est le quatrième moment centré et standardisé.

$m_4= \int_{-\infty}^\infty (\frac{z-\mu}{\sigma})^4P(z)dz$

Les distributions de Laplace, normale et uniforme montrées ont toutes une moyenne de 0 et une variance de 1, mais leur excès de kurtosis vaut 3, 0 et -1,2.

Les fonctions Python pour calculer les moments pourraient être :

# Import modules
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy
import scipy.stats as st
import mlgeo_synth

# One random generator for the whole notebook, seeded for reproducibility.
rng = np.random.default_rng(42)

def raw_moment(X, k, c=0):
    return ((X - c)**k).mean()


def central_moment(X, k):
    return raw_moment(X=X, k=k, c=X.mean())

2. Jeux de données géologiques [niveau 1]¶

Nous explorons la composition du granite en silice et en magnésium. Les données ont été collectées dans la base EarthChem.

# Load .csv data into a pandas dataframe
url = 'https://raw.githubusercontent.com/UW-MLGEO/MLGeo-dataset/main/data/EarthRocGranites.csv'
df = pd.read_csv(url)
df.head()

Le prétraitement des données est souvent nécessaire et, surtout, il est critique de consigner chaque étape de traitement appliquée aux données brutes. Ne modifiez pas le fichier de données original ; consignez plutôt les étapes de traitement. Ci-dessous, nous supprimons les lignes contenant des NaN (not a number).

df = df.dropna()    # remove rows with NaN values
df.head()

Pandas inclut des méthodes pour rapporter les statistiques de base des données. Utilisez la méthode describe du DataFrame.

df.describe()

# Now, let's visualize the histograms of silica and magnesium

# Create a subplot with two histograms side by side
fig, axes = plt.subplots(1, 2, figsize=(10, 4))  # 1 row, 2 columns

# Plot the histograms for each column
axes[0].hist(df['SIO2(WT%)'], bins=60, color='black')
axes[0].set_xlabel('SiO$_2$, wt%')
axes[0].set_ylabel('Count')
axes[0].set_xlim([40, 100])

axes[1].hist(df['MGO(WT%)'], bins=100, color='black')
axes[1].set_xlabel('MgO, wt%')
axes[1].set_ylabel('Count')
# Note these xlims -> the data largely [but not completely!] sit between 0 and 10 wt%
axes[1].set_xlim([0, 10])

# Add spacing between subplots
plt.tight_layout()
plt.show()

# One more plot: a scatter of SiO2 vs. MgO
plt.scatter(df['SIO2(WT%)'], df['MGO(WT%)'], c='red', alpha=0.125)
ax = plt.gca()
ax.set_xlim([0, 100])
ax.set_xlabel('SiO$_2$, wt%')
ax.set_ylim([0, 100])
ax.set_ylabel('MgO, wt%')
ax.set_aspect('equal')

Calculons maintenant les moments pour SiO₂ avec les fonctions que nous avons définies plus haut.

# The mean:
print(f'The mean is: {raw_moment(df["SIO2(WT%)"], 1):4.2f}')

# Variance:
print(f'The variance is: {central_moment(df["SIO2(WT%)"], 2):4.2f}')

# Skewness:
skewness = central_moment(df["SIO2(WT%)"], 3) / central_moment(df["SIO2(WT%)"], 2) ** (3/2)
print(f'The skewness is: {skewness:4.2f}')

# Kurtosis:
kurtosis_value = central_moment(df['SIO2(WT%)'], 4) / central_moment(df['SIO2(WT%)'], 2) ** 2
print(f'The kurtosis is: {kurtosis_value:4.2f}')

The mean is: 72.11
The variance is: 16.84
The skewness is: -1.75
The kurtosis is: 13.67

# We can also just use pandas (or numpy or scipy):
print('The mean is: %4.2f, the variance is: %4.2f, the skewness is: %4.2f, and the kurtosis is: %4.2f'
      % (df['SIO2(WT%)'].mean(), df['SIO2(WT%)'].var(), df['SIO2(WT%)'].skew(), df['SIO2(WT%)'].kurtosis()))

The mean is: 72.11, the variance is: 16.84, the skewness is: -1.75, and the kurtosis is: 10.67

Notez que pandas rapporte l’excès de kurtosis (distribution normale = 0), tandis que notre version central_moment rapporte le kurtosis brut (distribution normale = 3). Gardez trace de la convention qu’utilise chaque bibliothèque.

3. Distributions géoscientifiques¶

Exemple 1 : échantillonner la distribution normale

Application : simuler des variations journalières de température en un lieu donné au cours du temps. C’est un substitut synthétique, pas des observations : nous tirons d’une distribution normale avec des paramètres plausibles pour Seattle (moyenne annuelle autour de 11,5 °C, écart-type autour de 6 °C). Les vraies données de température ont une structure saisonnière qu’une seule distribution normale ne capture pas.

# Synthetic stand-in for daily mean temperature in Seattle.
# These are NOT observations; the parameters are plausible values for Seattle.
mean_temp = 11.5  # annual mean temperature (C), plausible for Seattle
std_temp = 6.0    # standard deviation (C), plausible for Seattle

temperatures = rng.normal(loc=mean_temp, scale=std_temp, size=1000)

plt.hist(temperatures, bins=30, color='skyblue', edgecolor='black')
plt.title('Synthetic Temperature Distribution (Normal, Seattle-like parameters)')
plt.xlabel('Temperature (C)')
plt.ylabel('Frequency')
plt.show()

Exemple 2 : magnitudes des séismes et loi de Gutenberg-Richter

Les magnitudes des séismes ne suivent pas une distribution log-normale. Elles suivent la loi de Gutenberg-Richter : le nombre de séismes $N$ de magnitude au moins $M$ vérifie

$\log_{10} N = a - bM$ ,

où $a$ fixe le taux global de sismicité et $b$ (la « valeur de b », b-value) le taux relatif des petits événements par rapport aux grands. Globalement, $b \approx 1$ : à chaque unité de magnitude en moins, il y a environ dix fois plus de séismes. La magnitude étant déjà une mesure logarithmique de la taille, cette loi signifie que les magnitudes des séismes suivent une distribution exponentielle au-dessus de la magnitude minimale du catalogue.

Nous générons un catalogue synthétique avec mlgeo_synth.gutenberg_richter_magnitudes, qui tire des magnitudes avec une valeur de b connue.

b_true = 1.0
m_min = 1.0
magnitudes = mlgeo_synth.gutenberg_richter_magnitudes(n=20000, b=b_true, m_min=m_min,
                                                      m_max=8.0, seed=42)
print(f'{len(magnitudes)} magnitudes between {magnitudes.min():.2f} and {magnitudes.max():.2f}')

20000 magnitudes between 1.00 and 5.64

Tracez la distribution fréquence-magnitude. Avec l’axe des comptes en échelle logarithmique, la loi de Gutenberg-Richter apparaît comme une droite de pente $-b$ : les comptes par classe de magnitude comme les comptes cumulés $N(\geq M)$ décroissent log-linéairement.

bins = np.arange(1.0, 8.1, 0.1)
counts, edges = np.histogram(magnitudes, bins=bins)
bin_centers = 0.5 * (edges[:-1] + edges[1:])

# Cumulative count of events with magnitude >= M
mags_sorted = np.sort(magnitudes)
n_cum = len(magnitudes) - np.arange(len(magnitudes))

fig, ax = plt.subplots(figsize=(7, 5))
ax.semilogy(bin_centers[counts > 0], counts[counts > 0], 'ks', ms=4,
            label='counts per 0.1 bin')
ax.semilogy(mags_sorted, n_cum, 'r-', lw=2, label=r'cumulative $N(\geq M)$')
ax.set_xlabel('Magnitude M')
ax.set_ylabel('Number of earthquakes (log scale)')
ax.set_title('Frequency-magnitude distribution (Gutenberg-Richter)')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
plt.show()

La droite sur l’axe des comptes en échelle logarithmique est la signature de la loi de Gutenberg-Richter.

Nous pouvons estimer la valeur de b à partir des données avec l’estimateur du maximum de vraisemblance (Aki, 1965) :

$\hat{b} = \dfrac{\log_{10}(e)}{\bar{M} - M_{min}}$ ,

où $\bar{M}$ est la magnitude moyenne du catalogue et $M_{min}$ la magnitude minimale de complétude.

b_est = np.log10(np.e) / (np.mean(magnitudes) - m_min)
print(f'Maximum-likelihood b-value estimate: {b_est:.3f} (true value: {b_true})')

Maximum-likelihood b-value estimate: 0.998 (true value: 1.0)

L’estimation retrouve la valeur de b que nous avons injectée. Sur les catalogues réels, le même estimateur fonctionne une fois que vous avez identifié la magnitude de complétude (la magnitude au-dessus de laquelle le réseau détecte tous les événements) ; en dessous, le catalogue manque les petits séismes et la droite s’infléchit.

Exemple 3 : distributions en loi puissance

Application en géosciences : les distributions en loi puissance s’observent dans les occurrences d’aléas naturels comme les glissements de terrain et les feux de forêt, où les petits événements sont fréquents mais les grands événements rares.

# Generating samples from a power-law distribution
a = 2.5  # Shape parameter (the larger, the steeper the fall-off)
size = 1000

power_law_data = (rng.pareto(a, size) + 1) * 10  # Shifted Pareto distribution

plt.hist(power_law_data, bins=50, color='red', edgecolor='black', log=True)
plt.title('Simulated Data from Power-Law Distribution')
plt.xlabel('Event Size')
plt.ylabel('Frequency (log scale)')
plt.show()

La distribution en loi puissance capture le comportement à queue lourde typique des processus géophysiques comme les glissements de terrain.

4. Exercice en classe¶

Choisissez l’une des distributions de cette leçon, tirez-en des échantillons, calculez les quatre premiers moments et comparez-les aux valeurs théoriques. Suivez les étapes de la cellule ci-dessous.

# In-class exercise: moments of a distribution.
#
# Step 1: pick a distribution and draw 5000 samples with the generator, e.g. one of:
#   samples = rng.normal(loc=11.5, scale=6.0, size=5000)
#   samples = mlgeo_synth.gutenberg_richter_magnitudes(n=5000, b=1.0, m_min=1.0, seed=42)
#   samples = (rng.pareto(2.5, 5000) + 1) * 10
#
# Step 2: compute the first four moments with the functions from Section 1:
#   mean     -> raw_moment(samples, 1)
#   variance -> central_moment(samples, 2)
#   skewness -> central_moment(samples, 3) / central_moment(samples, 2)**(3/2)
#   kurtosis -> central_moment(samples, 4) / central_moment(samples, 2)**2
#
# Step 3: compare with the theoretical values.
#   Normal(mu, sigma): mean mu, variance sigma^2, skewness 0, kurtosis 3.
#   Exponential (Gutenberg-Richter above m_min, rate beta = b*ln(10)):
#     mean m_min + 1/beta, variance 1/beta^2, skewness 2, kurtosis 9.
#
# Step 4: cross-check with scipy: st.skew(samples), st.kurtosis(samples, fisher=False).
#
# Step 5: repeat with only 100 samples. Which moments are most sensitive
#         to sample size, and why?