The Gumbel Distribution, Extreme Rainfall, and the Euler–Mascheroni Constant

"""
Histogram of Annual Maxima with Gumbel Fit (Method of Moments)

Overview
--------
This script simulates daily rainfall amounts for many years, computes the
annual maximum rainfall, and fits a Gumbel distribution to those maxima using
the method of moments (MoM). It then overlays the fitted Gumbel PDF on the
histogram of annual maxima and saves the figure as a PNG.

Why Gumbel?
-----------
The Gumbel distribution often appears as a limiting distribution for maxima
(e.g., yearly maximum daily rainfall), under broad conditions (Extreme Value
Theory).

Model Assumptions
-----------------
- Daily rainfall is modeled as Exponential(mean = 10 mm).
  This captures "many small, few large" rainfall days.
- 365 days per year.
- We simulate `YEARS` independent years to get a smooth histogram.

Gumbel(μ, β) Method-of-Moments
-------------------------------
For a random variable X ~ Gumbel(μ, β),
- E[X] = μ + γ β,  where γ is the Euler–Mascheroni constant
- Var[X] = (π^2 / 6) β^2
Thus, given sample mean m and sample variance s^2, the MoM estimators are:
- β̂ = sqrt(6 s^2) / π
- μ̂ = m − γ β̂

Outputs
-------
- PNG figure: 'gumbel_annual_max_hist_fit.png'
- The console prints fitted parameters μ̂ and β̂

Usage
-----
Run this file with Python 3. It requires only NumPy and Matplotlib.
"""

import numpy as np
import matplotlib.pyplot as plt

# ------------------- Configuration -------------------
SEED = 20251011
DAYS_PER_YEAR = 365
YEARS = 20000          # large for a smooth fit
DAILY_MEAN_MM = 10.0   # mean of exponential daily rainfall

# ------------------- Simulation ----------------------
rng = np.random.default_rng(SEED)

# Simulate daily rainfall: shape (YEARS, DAYS_PER_YEAR)
daily_rain = rng.exponential(DAILY_MEAN_MM, size=(YEARS, DAYS_PER_YEAR))

# Annual maxima: one maximum per simulated year
annual_max = daily_rain.max(axis=1)

# ------------------- Gumbel MoM fit ------------------
GAMMA = 0.5772156649015329

xbar = annual_max.mean()
s2 = annual_max.var(ddof=1)              # unbiased sample variance
beta_hat = np.sqrt(6.0 * s2) / np.pi     # β̂
mu_hat = xbar - GAMMA * beta_hat         # μ̂

print(f"Fitted Gumbel parameters (MoM): mu = {mu_hat:.6f}, beta = {beta_hat:.6f}")

# ------------------- Overlay PDF ---------------------
xs = np.linspace(annual_max.min(), annual_max.max(), 600)
z = (xs - mu_hat) / beta_hat
gumbel_pdf = (1.0 / beta_hat) * np.exp(-(z + np.exp(-z)))

# ------------------- Plot & Save ---------------------
plt.figure(figsize=(7,5))
plt.hist(annual_max, bins=60, density=True, alpha=0.6, label="Annual maxima")
plt.plot(xs, gumbel_pdf, label=f"Gumbel fit (μ={mu_hat:.2f}, β={beta_hat:.2f})")
plt.title("Annual Maximum Daily Rainfall\nHistogram with Gumbel PDF (method-of-moments fit)")
plt.xlabel("Max rainfall in a year (mm)")
plt.ylabel("Density")
plt.legend()

OUTPUT = "gumbel_annual_max_hist_fit.png"
plt.savefig(OUTPUT, dpi=200, bbox_inches="tight")
print(f"Saved figure to: {OUTPUT}")
# (Optional) plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path

# ===== USER SETTINGS =====
STATION_FILE = "RR_STAID012345.txt"  # change to your station file
START_YEAR = 1941                    # analysis start
MIN_DAYS_PER_YEAR = 330              # completeness (drops partial years like 2025)
OUT_DIR = Path("gumbel_out")
# =========================

GAMMA = 0.5772156649015329  # Euler–Mascheroni
OUT_DIR.mkdir(exist_ok=True, parents=True)

# ---- helpers ----
def detect_header_row(path: Path) -> int:
    lines = path.read_text(encoding="utf-8", errors="ignore").splitlines()
    for i, ln in enumerate(lines[:120]):
        if "SOUID" in ln and "DATE" in ln and "RR" in ln:
            return i
    return 0

def read_ecad_station(path: Path) -> pd.DataFrame:
    hdr = detect_header_row(path)
    df = pd.read_csv(path, skiprows=hdr, header=0)
    df.columns = [c.strip() for c in df.columns]
    keep = [c for c in df.columns if c in ("SOUID","DATE","RR","Q_RR")]
    df = df[keep]
    df["DATE"] = pd.to_datetime(df["DATE"].astype(str), format="%Y%m%d", errors="coerce")
    df["RR"]   = pd.to_numeric(df["RR"], errors="coerce")
    # RR is in 0.1 mm; -9999 means missing
    df["RR_mm"] = df["RR"].where(df["RR"] != -9999, np.nan) / 10.0
    # keep good-quality only (if flag present)
    if "Q_RR" in df.columns:
        df["Q_RR"] = pd.to_numeric(df["Q_RR"], errors="coerce")
        df = df[df["Q_RR"] == 0]
    df["year"] = df["DATE"].dt.year
    return df.dropna(subset=["DATE"])

def latest_complete_year(df: pd.DataFrame, min_days=MIN_DAYS_PER_YEAR):
    counts = df.groupby("year")["RR_mm"].apply(lambda s: s.notna().sum())
    complete = counts[counts >= min_days]
    return None if complete.empty else int(complete.index.max())

def gumbel_mom(sample):
    sample = np.asarray(sample, float)
    m  = float(sample.mean())
    s2 = float(sample.var(ddof=1))
    beta = np.sqrt(6.0 * s2) / np.pi
    mu   = m - GAMMA * beta
    return mu, beta

def gumbel_pdf(z, mu, beta):
    y = (z - mu) / beta
    return (1.0/beta) * np.exp(-(y + np.exp(-y)))

# ---- 1) Read & keep only complete years ----
df = read_ecad_station(Path(STATION_FILE))
if df.empty:
    raise SystemExit("No usable rows after filtering. Check the file.")

last_full = latest_complete_year(df, MIN_DAYS_PER_YEAR)
if last_full is None:
    raise SystemExit("No complete years under the chosen threshold.")

df = df[(df["year"] >= START_YEAR) & (df["year"] <= last_full)].copy()
if df.empty:
    raise SystemExit("No data left in 1941..latest complete year after filtering.")

# ---- 2) Annual maxima ----
annual_max = df.groupby("year")["RR_mm"].max(min_count=1).dropna()
years = annual_max.index.values
x = annual_max.values
t = years - years.min()
n = len(x)

# ---- 3) Fit Gumbel (no trend) + compute OLS trend for series plot ----
mu0, beta0 = gumbel_mom(x)
A = np.vstack([np.ones_like(t), t]).T
coef, *_ = np.linalg.lstsq(A, x, rcond=None)
a_hat, b_hat = float(coef[0]), float(coef[1])

# ---- 4) PLOTS (only two) ----
yr_label = f"{years.min()}–{years.max()} (complete years)"

# 4.1 Histogram + Gumbel fit (no trend)
xs = np.linspace(x.min() - 5, x.max() + 5, 400)
plt.figure(figsize=(7,5))
plt.hist(x, bins=20, density=True, alpha=0.6, edgecolor="black", label="Annual maxima")
plt.plot(xs, gumbel_pdf(xs, mu0, beta0), linewidth=2, label="Gumbel fit (no trend)")
plt.xlabel("Annual maximum daily precipitation (mm)")
plt.ylabel("Density")
plt.title(f"Histogram with Gumbel fit\n{yr_label}")
plt.legend()
plt.tight_layout()
plt.savefig(OUT_DIR / "hist_gumbel_notrend.png", dpi=160)
plt.close()

# 4.2 Time series with OLS trend
plt.figure(figsize=(9,4))
plt.plot(years, x, marker="o", linestyle="-", label="Annual max (Rx1day)")
plt.plot(years, a_hat + b_hat * t, linewidth=2, label=f"OLS trend (slope {b_hat:.2f} mm/yr)")
plt.xlabel("Year")
plt.ylabel("Annual maximum daily precipitation (mm)")
plt.title(f"Annual maxima with linear trend\n{yr_label}")
plt.legend()
plt.tight_layout()
plt.savefig(OUT_DIR / "series_trend.png", dpi=160)
plt.close()

print("Saved:", OUT_DIR / "hist_gumbel_notrend.png", "and", OUT_DIR / "series_trend.png")
print(f"Gumbel (no trend): mu={mu0:.2f} mm, beta={beta0:.2f} mm; OLS slope={b_hat:.3f} mm/yr")