Add ensemble calibration diagnostics to inst/ilamb

Adds a module of ensemble calibration diagnostics (rank histogram,
spread-skill ratio, coverage, reliability) that assess whether an
ensemble's spread is well calibrated against observations, complementing
the existing ensemble-mean benchmarking. Includes synthetic-ensemble
tests, an example driver, and an optional plotting script.

Author: Tejas7007

CHANGELOG.md

@@ -9,6 +9,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 ## Unreleased
 
 ### Added
+- Added ensemble calibration diagnostics to `inst/ilamb/` in PEcAn.benchmark (rank histogram, spread-skill ratio, coverage, reliability) for assessing whether an ensemble's spread is well calibrated against observations, complementing the ensemble-mean benchmarking.
 - Added `make_scorecard.sh` and documentation to `inst/ilamb/` in PEcAn.benchmark for generating and serving the ILAMB HTML scorecard.
 - New function `PEcAn.utils::netcdf2df()` flattens all dims and vars of a netCDF into a dataframe,
     with units attached as an attribute.

modules/benchmark/inst/ilamb/calibration/README.md

@@ -0,0 +1,100 @@
+# Ensemble Calibration Diagnostics
+
+Diagnostics for assessing whether an ensemble's *spread* is trustworthy, not
+just whether its *mean* is accurate. This complements the ILAMB benchmarking
+already in this directory, which scores ensemble-mean skill. Here we ask a
+different question: when a model reports uncertainty through an ensemble, is that
+uncertainty well calibrated against observations?
+
+## Why this matters
+
+Carbon reanalyses and multi-model ensembles increasingly report ensemble spread
+as an uncertainty estimate. A skillful mean is not enough. If the spread is too
+narrow the product is overconfident and the observation will routinely fall
+outside the ensemble range. If it is too wide the product is underconfident.
+These diagnostics quantify where an ensemble sits on that spectrum.
+
+## The diagnostics
+
+All functions take `members` (a NumPy array of shape `[n_members, ...]`) and
+`obs` (an array matching the per-cell shape). Missing values in either are
+handled by masking to cells where the observation and all members are present.
+
+- **`rank_histogram(members, obs)`** — the histogram of where each observation
+  falls within its sorted ensemble. A flat histogram indicates calibration, a
+  U shape indicates overconfidence (the observation keeps landing outside the
+  ensemble), and a dome indicates underconfidence.
+
+- **`spread_skill(members, obs)`** — the ratio of ensemble spread to ensemble
+  error, with a small-sample correction. A ratio near 1 indicates calibration,
+  well below 1 indicates overconfidence, and above 1 indicates underconfidence.
+  This ratio is the metric to use when comparing ensembles of different sizes,
+  because it carries the sqrt((n+1)/n) correction (see the note below).
+
+- **`coverage(members, obs, interval=None)`** — the fraction of observations
+  that fall inside the ensemble range, or inside a central interval if one is
+  given (for example `interval=0.9` for the central 90 percent). For a
+  calibrated ensemble this should match the expected value, which for the full
+  range is (n-1)/(n+1).
+
+- **`reliability(members, obs, threshold)`** — for a chosen threshold, compares
+  the ensemble's forecast probability of exceeding it against the observed
+  frequency. A calibrated ensemble lies on the diagonal.
+
+## A note on comparing ensembles of different sizes
+
+Fixed-interval coverage is not directly comparable across ensembles with very
+different member counts. A small ensemble cannot resolve, say, a 90 percent
+central interval, so its coverage is legitimately below the nominal target even
+when the ensemble is calibrated. When comparing a large ensemble against a small
+one, prefer the bias-corrected `spread_skill` ratio, which is designed to be
+fair across sizes. The test suite documents this behavior explicitly.
+
+## Usage
+
+```python
+import numpy as np
+import ensemble_calibration as ec
+
+# members: shape [n_members, n_cells], obs: shape [n_cells]
+ratio = ec.spread_skill(members, obs)["ratio"]
+cov   = ec.coverage(members, obs, interval=0.9)["coverage"]
+counts, edges = ec.rank_histogram(members, obs)
+
+if ratio < 0.8:
+    print("ensemble is overconfident relative to these observations")
+elif ratio > 1.2:
+    print("ensemble is underconfident")
+else:
+    print("ensemble spread is well sized")
+```
+
+See `example_calibration.py` for a full example that loads model output and
+benchmark fields and reports the diagnostics for several variables.
+
+## Testing
+
+`test_ensemble_calibration.py` validates every diagnostic against synthetic
+ensembles whose calibration is known by construction. An ensemble built to be
+overconfident must be flagged as overconfident, a calibrated one must read as
+calibrated, and so on. The suite also checks that the calibration verdict is
+stable across member counts from 8 to 100, and that the diagnostics handle
+ragged, partially missing model output without error. Run with:
+
+```
+python test_ensemble_calibration.py
+```
+
+## Visualizing the results
+
+`plot_calibration.py` produces a two-panel figure summarizing the calibration
+comparison: panel (a) the bias-corrected spread-to-error ratio per variable,
+and panel (b) the fraction of observations inside the ensemble 90 percent band.
+Both panels are comparable across ensembles of different sizes. An example of
+the output is shown below.
+
+![Calibration comparison](calibration_figure.png)
+
+The plotting script additionally requires matplotlib and xarray. The core
+diagnostics in `ensemble_calibration.py` remain pure NumPy and have no plotting
+or file-format dependencies.

modules/benchmark/inst/ilamb/calibration/calibration_figure.png

@@ -0,0 +1,176 @@
+"""
+Ensemble calibration diagnostics for PEcAn benchmarking.
+
+The multi-model comparison in the ILAMB benchmarking evaluates ensemble MEANS:
+is the average field accurate? These diagnostics ask a different question, is
+the ensemble SPREAD the right size? An ensemble can have an accurate mean yet
+be overconfident (spread too narrow for its error) or underconfident (spread
+too wide).
+
+Four diagnostics are provided:
+
+  rank_histogram   For each grid cell, the rank of the observation among the
+                   sorted ensemble members. Aggregated over cells, a flat
+                   histogram indicates calibration; a U shape indicates the
+                   ensemble is overconfident (truth falls outside the members
+                   too often); a dome indicates underconfidence.
+
+  spread_skill     Ensemble spread (standard deviation across members) versus
+                   ensemble error (absolute deviation of the mean from the
+                   observation). For a well-calibrated ensemble these are
+                   comparable in magnitude; spread much smaller than error
+                   indicates overconfidence.
+
+  coverage         Fraction of grid cells where the observation falls within
+                   the ensemble range (min to max, or a central interval).
+                   Reported against the value expected for a calibrated
+                   ensemble of the same size.
+
+  reliability      For an event defined by a threshold, the ensemble forecast
+                   probability (fraction of members exceeding it) binned
+                   against the observed frequency of the event. A calibrated
+                   ensemble lies on the diagonal.
+
+A stochastic-member ensemble (such as a PEcAn data-assimilation reanalysis) is
+a true ensemble, so these diagnostics apply to it directly. A multi-model
+ensemble (different models, not draws from one distribution) can be scored the
+same way, but the diagnostics then test whether that inter-model spread behaves
+like a calibrated uncertainty distribution, which is a distinct and largely
+untested question.
+
+All functions take numpy arrays. Members are stacked on the first axis:
+`members` has shape (n_members, ...) and `obs` has shape (...), matching the
+trailing dimensions. NaNs (missing cells) are handled by masking to cells
+where the observation and all members are finite.
+"""
+
+import numpy as np
+
+
+def _valid_mask(members, obs):
+    """Cells where obs and every member are finite."""
+    return np.isfinite(obs) & np.all(np.isfinite(members), axis=0)
+
+
+def rank_histogram(members, obs, n_bins=None):
+    """
+    Rank of the observation within the ensemble, per valid cell.
+
+    Returns (counts, edges): `counts` has length n_members + 1 (the number of
+    possible ranks), giving how many observations fell at each rank. Ties are
+    broken by adding a small uniform jitter, the standard treatment so ties do
+    not artificially pile up at one rank. The jitter uses a fixed seed so the
+    result is reproducible.
+    """
+    members = np.asarray(members, dtype=float)
+    obs = np.asarray(obs, dtype=float)
+    n = members.shape[0]
+    mask = _valid_mask(members, obs)
+
+    m = members[:, mask]
+    o = obs[mask]
+
+    rng = np.random.default_rng(0)
+    jitter = rng.uniform(-1e-9, 1e-9, size=m.shape)
+    mj = m + jitter
+    oj = o + rng.uniform(-1e-9, 1e-9, size=o.shape)
+    ranks = (mj < oj).sum(axis=0)
+
+    counts = np.bincount(ranks, minlength=n + 1)
+    edges = np.arange(n + 2)
+    return counts, edges
+
+
+def spread_skill(members, obs):
+    """
+    Ensemble spread versus error.
+
+    Returns a dict with:
+      spread     mean over cells of the ensemble standard deviation
+      rmse       root mean square error of the ensemble mean vs obs
+      ratio      spread / rmse (near 1 for a well-calibrated ensemble;
+                 much less than 1 indicates overconfidence)
+    The spread is scaled by sqrt((n+1)/n) so it is comparable to the error of
+    the mean for a finite ensemble.
+    """
+    members = np.asarray(members, dtype=float)
+    obs = np.asarray(obs, dtype=float)
+    n = members.shape[0]
+    mask = _valid_mask(members, obs)
+
+    m = members[:, mask]
+    o = obs[mask]
+    mean = m.mean(axis=0)
+    std = m.std(axis=0, ddof=1)
+
+    spread = float(np.sqrt((n + 1) / n) * std.mean())
+    rmse = float(np.sqrt(np.mean((mean - o) ** 2)))
+    ratio = spread / rmse if rmse > 0 else float("nan")
+    return {"spread": spread, "rmse": rmse, "ratio": ratio}
+
+
+def coverage(members, obs, interval=None):
+    """
+    Fraction of cells where obs falls within the ensemble range.
+
+    interval=None uses the full member min-max. Otherwise pass a central
+    fraction (e.g. 0.9) to use that central interval via member quantiles.
+    Returns a dict with the observed coverage and, for the full range, the
+    value expected for a calibrated ensemble of this size: (n-1)/(n+1).
+    """
+    members = np.asarray(members, dtype=float)
+    obs = np.asarray(obs, dtype=float)
+    n = members.shape[0]
+    mask = _valid_mask(members, obs)
+
+    m = members[:, mask]
+    o = obs[mask]
+
+    if interval is None:
+        lo = m.min(axis=0)
+        hi = m.max(axis=0)
+        expected = (n - 1) / (n + 1)
+    else:
+        q = (1 - interval) / 2
+        lo = np.quantile(m, q, axis=0)
+        hi = np.quantile(m, 1 - q, axis=0)
+        expected = interval
+
+    inside = float(np.mean((o >= lo) & (o <= hi)))
+    return {"coverage": inside, "expected": expected, "n_members": n}
+
+
+def reliability(members, obs, threshold, n_bins=10):
+    """
+    Reliability of the ensemble probability for an event obs > threshold.
+
+    The forecast probability at each cell is the fraction of members exceeding
+    the threshold. Cells are grouped into probability bins; for each bin the
+    mean forecast probability and the observed event frequency are returned.
+    A calibrated ensemble lies on the diagonal (forecast == observed).
+
+    Returns a dict with bin_prob (mean forecast probability per bin),
+    obs_freq (observed frequency per bin), and count (cells per bin).
+    """
+    members = np.asarray(members, dtype=float)
+    obs = np.asarray(obs, dtype=float)
+    mask = _valid_mask(members, obs)
+
+    m = members[:, mask]
+    o = obs[mask]
+    p = (m > threshold).mean(axis=0)
+    y = (o > threshold).astype(float)
+
+    edges = np.linspace(0, 1, n_bins + 1)
+    idx = np.clip(np.digitize(p, edges) - 1, 0, n_bins - 1)
+
+    bin_prob = np.full(n_bins, np.nan)
+    obs_freq = np.full(n_bins, np.nan)
+    count = np.zeros(n_bins, dtype=int)
+    for b in range(n_bins):
+        sel = idx == b
+        count[b] = sel.sum()
+        if count[b] > 0:
+            bin_prob[b] = p[sel].mean()
+            obs_freq[b] = y[sel].mean()
+    return {"bin_prob": bin_prob, "obs_freq": obs_freq, "count": count}