import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
AMBER="#d97706"; TEAL="#0d9488"; INK="#1a2138"; GRID="#e6e9f2"; PINK="#db2777"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try:    df = pd.read_csv("../../data/meteorology_regional_rainfall.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"meteorology_regional_rainfall.csv")
print("loaded:", df.shape)
df.head()

loaded: (1000, 4)

x = df["rainfall_inches"]
m, v = x.mean(), x.var()
alpha = m**2 / v          # shape
beta  = v / m             # scale
print(f"observations : {len(x):,}")
print(f"mean {m:.3f}, variance {v:.3f}, skewness {stats.skew(x):.3f} (positive -> right-skewed)")
print(f"method-of-moments shape alpha = mean^2/var = {alpha:.3f}")
print(f"method-of-moments scale beta  = var/mean   = {beta:.3f}")
print(f"check: mean alpha*beta = {alpha*beta:.3f}, var alpha*beta^2 = {alpha*beta**2:.3f}")

observations : 1,000
mean 4.040, variance 5.573, skewness 1.202 (positive -> right-skewed)
method-of-moments shape alpha = mean^2/var = 2.929
method-of-moments scale beta  = var/mean   = 1.379
check: mean alpha*beta = 4.040, var alpha*beta^2 = 5.573

fig,ax=plt.subplots(figsize=(7,3.6))
ax.hist(x,bins=45,density=True,color=AMBER,alpha=0.7,label="observed rainfall")
xs=np.linspace(0,x.max(),300); ax.plot(xs,stats.gamma.pdf(xs,alpha,scale=beta),color=TEAL,lw=2.5,label=f"Gamma(shape={alpha:.2f}, scale={beta:.2f})")
ax.axvline(m,color=PINK,ls="--",lw=2,label=f"mean = {m:.1f} in")
ax.set_xlabel("seasonal rainfall (inches)"); ax.set_ylabel("density"); ax.set_title("Rainfall: observed vs gamma"); ax.legend()
plt.tight_layout(); plt.show()

# Formal goodness-of-fit for the Gamma: KS test + a gamma QQ plot.
ksg = stats.kstest(x, "gamma", args=(alpha, 0, beta))
print(f"Kolmogorov-Smirnov vs Gamma: D = {ksg.statistic:.3f}, p = {ksg.pvalue:.3f}  -> "
      + ("cannot reject the fit (Gamma is a good model)" if ksg.pvalue>0.05 else "reject: not gamma"))
fig, ax = plt.subplots(figsize=(4.8,3.8))
stats.probplot(x, dist=stats.gamma, sparams=(alpha, 0, beta), plot=ax)
ax.get_lines()[0].set_color(AMBER); ax.get_lines()[0].set_markersize(3)
ax.set_title("Gamma QQ plot (points on the line = good fit)"); plt.tight_layout(); plt.show()

Kolmogorov-Smirnov vs Gamma: D = 0.019, p = 0.869  -> cannot reject the fit (Gamma is a good model)

g = stats.gamma(alpha, scale=beta)
for thr in [8,10,12,15]:
    p_exceed = g.sf(thr)
    rp = 1/p_exceed if p_exceed>0 else float("inf")
    print(f"P(rain > {thr:2d} in) = {p_exceed:.4f}   return period ~ 1-in-{rp:.0f} seasons  (observed {(x>thr).mean():.4f})")
print(f"\n95th percentile (a wet season) = {g.ppf(0.95):.2f} inches")
print(f"99th percentile (a flood year)  = {g.ppf(0.99):.2f} inches")

P(rain >  8 in) = 0.0664   return period ~ 1-in-15 seasons  (observed 0.0690)
P(rain > 10 in) = 0.0225   return period ~ 1-in-44 seasons  (observed 0.0210)
P(rain > 12 in) = 0.0072   return period ~ 1-in-139 seasons  (observed 0.0060)
P(rain > 15 in) = 0.0012   return period ~ 1-in-832 seasons  (observed 0.0010)

95th percentile (a wet season) = 8.54 inches
99th percentile (a flood year)  = 11.43 inches

rng = np.random.default_rng(56)
# sum of 3 exponentials (alpha~3) vs the fitted gamma
sum3 = rng.exponential(beta, size=(300_000,3)).sum(axis=1)
print(f"sum of 3 Exponential(scale={beta:.2f}): mean {sum3.mean():.3f}, var {sum3.var():.3f}")
print(f"Gamma(shape=3, scale={beta:.2f}):       mean {stats.gamma.mean(3,scale=beta):.3f}, var {stats.gamma.var(3,scale=beta):.3f}")
print(f"fitted gamma shape = {alpha:.2f}  (near 3 -> ~3 stacked rainfall pulses)")

sum of 3 Exponential(scale=1.38): mean 4.140, var 5.663
Gamma(shape=3, scale=1.38):       mean 4.138, var 5.709
fitted gamma shape = 2.93  (near 3 -> ~3 stacked rainfall pulses)

	station_id	year	season	rainfall_inches
0	STN_MIDWEST	1961	Spring	12.102
1	STN_MIDWEST	1929	Autumn	1.555
2	STN_ATLANTIC	2006	Autumn	6.018
3	STN_ATLANTIC	1946	Summer	1.454
4	STN_ATLANTIC	1936	Winter	2.016

⚙️ Setup & data¶