"""The verdict for "large earthquakes are becoming more frequent".

Committed 2026-09-11 (session 36) before any yearly count was fetched. The
data this function will be given is `yearly.csv` at minmag 7.0, 1973-2025.

Model: a quasi-Poisson regression of the annual count on the year. Earthquake
counts are overdispersed (aftershock sequences put several M7s in one year), so
the Poisson standard error is scaled by the square root of the Pearson
dispersion, floored at 1 so that an under-dispersed sample cannot make the
interval narrower than Poisson.

Verdict, on the change in the annual rate per decade and its 95 % interval:
  survived      the interval's lower bound is above zero (the rate rose);
  killed        otherwise, if the upper bound is below +10 % a decade
                (no rise, and a rise of the size a person would notice is
                ruled out: +10 % a decade is +66 % over the 53 years);
  inconclusive  otherwise (no rise detected, and a meaningful one not excluded).
"""

import numpy as np
import statsmodels.api as sm

MEANINGFUL = 0.10  # per decade


def trend(years, counts):
    years = np.asarray(years, dtype=float)
    counts = np.asarray(counts, dtype=float)
    x = (years - 2000.0) / 10.0  # decades, so the slope is per decade
    fit = sm.GLM(counts, sm.add_constant(x), family=sm.families.Poisson()).fit()
    mu = fit.fittedvalues
    phi = max(1.0, float(np.sum((counts - mu) ** 2 / mu) / (len(counts) - 2)))
    b = float(fit.params[1])
    se = float(fit.bse[1]) * np.sqrt(phi)
    return {
        "change": float(np.expm1(b)),
        "lo": float(np.expm1(b - 1.96 * se)),
        "hi": float(np.expm1(b + 1.96 * se)),
        "phi": phi,
        "n_years": len(counts),
        "mean_per_year": float(counts.mean()),
    }


def verdict(t):
    if t["lo"] > 0:
        return "survived"
    if t["hi"] < MEANINGFUL:
        return "killed"
    return "inconclusive"