2.7 Consideraciones estadísticas para los datos geocientíficos y el ruido

Distribuciones canónicas en los datos geocientíficos¶

Distribución normal: se usa para variables como la temperatura, las variaciones del nivel del mar o la velocidad del viento, cuyos valores se distribuyen simétricamente alrededor de una media.
Distribución log-normal: se observa en fenómenos cuyos valores no pueden ser negativos y muestran colas largas, como la intensidad de la lluvia, el caudal de los ríos, el tamaño de grano y la permeabilidad.
Distribución exponencial: se aplica al modelar los intervalos de tiempo entre eventos, como el tiempo entre sismos. Las magnitudes sísmicas siguen una forma exponencial relacionada, la ley de Gutenberg-Richter, que trabajamos más abajo.
Distribución de ley de potencias: aparece en eventos raros y de gran escala, como deslizamientos de tierra e incendios forestales, donde los eventos pequeños son comunes y los grandes son raros.

🖥️ Diapositivas — Sesión 07 (mié 14 oct)

1. Características estadísticas¶

Sea $P(z)$ la distribución de los datos $z$ .

La media¶

Imagen tomada de este blog.

La media es la suma de los valores dividida por el número de puntos de datos. Es el primer momento crudo de una distribución. $\mu = \int_{-\infty}^\infty zP(z)dz$ , donde $z$ es el valor de los datos (el intervalo o bin) y $P(z)$ es la distribución de los datos.

La varianza¶

La varianza es el segundo momento centralizado. Centralizado significa que la distribución se desplaza alrededor de la media. Calcula qué tan dispersa es una distribución.

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

La desviación estándar es la raíz cuadrada de la varianza, σ. Una varianza alta indica una distribución ancha.

La asimetría¶

La asimetría (skewness) es el tercer momento estandarizado. El momento estandarizado se escala por la desviación estándar. Mide el tamaño relativo de las dos colas de la distribución.

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

Con el exponente cúbico, es posible que la asimetría sea negativa.

Imagen tomada de este blog.

Una distribución con asimetría positiva es aquella donde la mayor parte del peso está al final de la distribución. Una distribución con asimetría negativa es aquella donde la mayor parte del peso está al principio de la distribución.

Curtosis¶

La curtosis mide el tamaño combinado de las dos colas en relación con la distribución completa. Es el cuarto momento centralizado y estandarizado.

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

kurtosis Las distribuciones de Laplace, normal y uniforme mostradas tienen todas media 0 y varianza 1, pero su exceso de curtosis es 3, 0 y -1.2.

Las funciones de Python para calcular los momentos podrían ser:

# 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. Conjuntos de datos geológicos [Nivel 1]¶

Exploramos la composición del granito en términos de su contenido de sílice y de magnesio. Los datos fueron recolectados de la base de datos 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()

El preprocesamiento de los datos suele ser necesario y, sobre todo, es crítico registrar cada paso de procesamiento aplicado a los datos crudos. No cambie el archivo de datos original; en su lugar, registre los pasos de procesamiento. Abajo eliminamos las filas con NaN (not a number, no es un número).

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

Pandas incluye métodos para reportar estadísticas básicas de los datos. Use el método describe del 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')

Ahora calculemos los momentos del SiO₂ con las funciones que definimos arriba.

# 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

Observe que pandas reporta el exceso de curtosis (distribución normal = 0), mientras que nuestra versión central_moment reporta la curtosis simple (distribución normal = 3). Lleve el registro de qué convención usa cada biblioteca.

3. Distribuciones geocientíficas¶

Ejemplo 1: muestreo de la distribución normal

Aplicación: simular las variaciones diarias de temperatura en una ubicación específica a lo largo del tiempo. Este es un sustituto sintético, no observaciones: extraemos de una distribución normal con parámetros plausibles para Seattle (media anual cercana a 11.5 C, desviación estándar cercana a 6 C). Los datos reales de temperatura tienen una estructura estacional que una sola distribución normal no captura.

# 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()

Ejemplo 2: magnitudes sísmicas y la ley de Gutenberg-Richter

Las magnitudes de los sismos no siguen una distribución log-normal. Siguen la ley de Gutenberg-Richter: el número de sismos $N$ con magnitud al menos $M$ satisface

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

donde $a$ fija la tasa general de sismicidad y $b$ (el «valor b») fija la tasa relativa de eventos pequeños frente a grandes. Globalmente, $b \approx 1$ : por cada unidad que baja la magnitud hay unas diez veces más sismos — es lo que se observa, por ejemplo, en los catálogos del Servicio Sismológico Nacional (SSN) de México. Como la magnitud ya es una medida logarítmica del tamaño, esta ley significa que las magnitudes sísmicas siguen una distribución exponencial por encima de la magnitud mínima del catálogo.

Generamos un catálogo sintético con mlgeo_synth.gutenberg_richter_magnitudes, que extrae magnitudes con un valor b conocido.

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

Grafique la distribución frecuencia-magnitud. Con el eje de conteos en escala logarítmica, la ley de Gutenberg-Richter aparece como una línea recta de pendiente $-b$ : tanto los conteos por intervalo de magnitud como los conteos acumulados $N(\geq M)$ decaen de forma log-lineal.

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 línea recta sobre el eje de conteos en escala logarítmica es la firma de la ley de Gutenberg-Richter.

Podemos estimar el valor b a partir de los datos con el estimador de máxima verosimilitud (Aki 1965):

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

donde $\bar{M}$ es la magnitud media del catálogo y $M_{min}$ es la magnitud mínima de completitud.

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)

La estimación recupera el valor b que pusimos. En catálogos reales, el mismo estimador funciona una vez que usted identificó la magnitud de completitud (la magnitud por encima de la cual la red detecta todos los eventos); por debajo de ella, el catálogo pierde los sismos pequeños y la línea se dobla.

Ejemplo 3: distribuciones de ley de potencias

Aplicación en las geociencias: las distribuciones de ley de potencias se observan en la ocurrencia de amenazas naturales como deslizamientos de tierra e incendios forestales, donde los eventos pequeños son comunes pero los grandes son raros.

# 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 distribución de ley de potencias captura el comportamiento de colas pesadas típico de procesos geofísicos como los deslizamientos de tierra.

4. Ejercicio en clase¶

Elija una de las distribuciones de esta lección, extraiga muestras de ella, calcule los primeros cuatro momentos y compárelos con los valores teóricos. Siga los pasos de la celda de abajo.

# 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?