Fitting Log-Pearson III Flood Frequency Curves in Python

import logging

import numpy as np
import pandas as pd
from scipy import stats

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


def frequency_factor(p_exceed: float, skew: float) -> float:
    """
    Pearson III frequency factor via the Wilson-Hilferty approximation.

    Reduces to the standard normal deviate at zero skew, which is the
    log-normal case, and is accurate across the skew range seen in practice.
    """
    z = stats.norm.ppf(1.0 - p_exceed)
    if abs(skew) < 1e-6:
        return float(z)
    k = skew / 6.0
    return float((2.0 / skew) * (((z - k) * k + 1.0) ** 3 - 1.0))


def station_skew_mse(station_skew: float, n: int) -> float:
    """
    Mean square error of a station skew estimate, as a function of record
    length. Short records give a large MSE, which is what pulls the weighted
    skew toward the regional value.
    """
    a = (-0.33 + 0.08 * abs(station_skew) if abs(station_skew) <= 0.9
         else -0.52 + 0.30 * abs(station_skew))
    b = (0.94 - 0.26 * abs(station_skew) if abs(station_skew) <= 1.5 else 0.55)
    return float(10 ** (a - b * np.log10(n / 10.0)))


def fit_lp3(
    peaks: pd.Series,
    regional_skew: float | None = None,
    regional_mse: float = 0.302,
    return_periods=(2, 5, 10, 25, 50, 100, 200, 500),
    confidence: float = 0.90,
) -> dict:
    """
    Fit Log-Pearson III and return quantiles, plotting positions and diagnostics.
    """
    s = pd.Series(peaks).dropna()
    s = s[s > 0]
    n = len(s)
    if n < 10:
        raise ValueError(f"{n} usable years is too short for a frequency fit")

    logs = np.log10(s.to_numpy(dtype=float))
    mean_log, std_log = float(logs.mean()), float(logs.std(ddof=1))
    g_station = float(stats.skew(logs, bias=False))

    skew, weight = g_station, 1.0
    if regional_skew is not None:
        mse_station = station_skew_mse(g_station, n)
        weight = regional_mse / (regional_mse + mse_station)
        skew = weight * g_station + (1.0 - weight) * regional_skew
        log.info("Station skew %.3f (MSE %.3f, weight %.2f), regional %.3f "
                 "→ weighted %.3f", g_station, mse_station, weight,
                 regional_skew, skew)

    z_conf = stats.norm.ppf(0.5 + confidence / 2.0)
    rows = []
    for T in return_periods:
        p = 1.0 / T
        kt = frequency_factor(p, skew)
        log_q = mean_log + kt * std_log
        # The quantile's standard error grows with |Kt| and shrinks with n,
        # which is why the band flares at long return periods on short records.
        se = std_log * np.sqrt((1.0 + 0.5 * kt ** 2) / n)
        rows.append({
            "return_period": T,
            "kt": kt,
            "estimate": 10 ** log_q,
            "lower": 10 ** (log_q - z_conf * se),
            "upper": 10 ** (log_q + z_conf * se),
        })

    quantiles = pd.DataFrame(rows).set_index("return_period")

    # --- Plotting positions for the visual check. Weibull is unbiased for the
    # exceedance probability, which is what we are plotting against. ---
    ranked = np.sort(s.to_numpy(dtype=float))[::-1]
    ranks = np.arange(1, n + 1)
    plotting = pd.DataFrame({
        "peak": ranked,
        "exceedance": ranks / (n + 1.0),
        "return_period": (n + 1.0) / ranks,
    })

    q100 = quantiles.loc[100]
    log.info("n=%d, mean(log)=%.4f, sd(log)=%.4f, skew=%.3f", n, mean_log, std_log, skew)
    log.info("Q100 = %.1f  [%.1f, %.1f]  — band spans a factor of %.2f",
             q100.estimate, q100.lower, q100.upper, q100.upper / q100.lower)
    if n < 100:
        log.warning("Record is %d years; the 100-year estimate is an "
                    "extrapolation of %.1f× the record length", n, 100 / n)

    # --- Leave-one-out on the largest peak: if the answer moves a lot, the
    # result is one flood's opinion and must be reported as such. ---
    without_max = fit_lp3_quantile(s.drop(s.idxmax()), 100, regional_skew, regional_mse)
    shift = 100.0 * (without_max - q100.estimate) / q100.estimate
    log.info("Dropping the largest peak moves Q100 by %+.1f %%", shift)

    return {
        "n": n, "mean_log": mean_log, "std_log": std_log,
        "station_skew": g_station, "weighted_skew": skew, "skew_weight": weight,
        "quantiles": quantiles, "plotting_positions": plotting,
        "q100_sensitivity_pct": shift,
    }


def fit_lp3_quantile(peaks: pd.Series, T: int, regional_skew=None,
                     regional_mse: float = 0.302) -> float:
    """A single quantile, used by the leave-one-out check."""
    logs = np.log10(pd.Series(peaks).dropna().pipe(lambda x: x[x > 0]).to_numpy(float))
    n = len(logs)
    g = float(stats.skew(logs, bias=False))
    if regional_skew is not None:
        w = regional_mse / (regional_mse + station_skew_mse(g, n))
        g = w * g + (1 - w) * regional_skew
    return float(10 ** (logs.mean() + frequency_factor(1.0 / T, g) * logs.std(ddof=1)))


# --- Example usage ---
# result = fit_lp3(natural_peaks, regional_skew=-0.05)
# print(result["quantiles"].round(1))