diff --git a/docs/api/index.md b/docs/api/index.md index 1b1b9179..85a7bc56 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -104,6 +104,15 @@ Use `ml4t.diagnostic.metrics` for reusable metric and feature-statistic helpers. - maximum_drawdown - hit_rate - compute_forward_returns + - brier_score + - brier_decomposition + - BrierDecomposition + - log_loss + - reliability_table + - expected_calibration_error + - max_calibration_error + - calibration_slope + - CalibrationSlope ## Cross-Validation diff --git a/docs/methods/calibration.md b/docs/methods/calibration.md new file mode 100644 index 00000000..9b14853f --- /dev/null +++ b/docs/methods/calibration.md @@ -0,0 +1,188 @@ +# Probability Calibration + +A probability forecast is calibrated when events forecast at probability `p` +happen with frequency `p`. The forecast can be a classifier's output, an +analyst's probability, or a prediction-market price read as an implied +probability. Calibration is separate from discrimination: a forecaster that +always predicts the base rate is calibrated but useless, and a model that ranks +events well can still be systematically overconfident. + +`ml4t.diagnostic.metrics` provides proper scores, the Brier decomposition, a +reliability table with confidence intervals, calibration-error summaries, and a +calibration slope. All functions take NumPy arrays, Polars or pandas Series, +or sequences of forecasts in `[0, 1]` and outcomes in `{0, 1}`. Missing values +raise an error instead of being dropped. + +## Example: prices with a favorite-longshot bias + +The synthetic market below quotes prices in cents. The true win probability is +`expit(1.4 * logit(price))`, so longshots win less often than priced and +favorites more often. Markets come in groups of four that resolve on one +event, so outcomes within a group are identical. + +```python +import numpy as np + +from ml4t.diagnostic.metrics import ( + brier_decomposition, + brier_score, + calibration_slope, + expected_calibration_error, + log_loss, + max_calibration_error, + reliability_table, +) + +rng = np.random.default_rng(0) +n_events, markets_per_event = 12_500, 4 +event_price = np.clip(np.round(rng.beta(0.6, 0.6, n_events), 2), 0.01, 0.99) +true_prob = 1 / (1 + np.exp(-1.4 * np.log(event_price / (1 - event_price)))) +event_outcome = (rng.uniform(size=n_events) < true_prob).astype(int) + +price = np.repeat(event_price, markets_per_event) +outcome = np.repeat(event_outcome, markets_per_event) +event_id = np.repeat(np.arange(n_events), markets_per_event) + +print(f"Brier score: {brier_score(price, outcome):.4f}") +print(f"Log loss: {log_loss(price, outcome):.4f}") + +decomposition = brier_decomposition(price, outcome) +print( + f"reliability={decomposition.reliability:.4f} " + f"resolution={decomposition.resolution:.4f} " + f"uncertainty={decomposition.uncertainty:.4f}" +) + +table = reliability_table(price, outcome, clusters=event_id, n_bootstrap=200, seed=0) +print(table.select("bin_lower", "bin_upper", "n", "n_clusters", + "mean_forecast", "observed_frequency", "ci_lower", "ci_upper")) + +print(f"ECE: {expected_calibration_error(price, outcome):.4f}") +print(f"MCE: {max_calibration_error(price, outcome):.4f}") + +slope = calibration_slope(price, outcome, spec="logit", clusters=event_id) +print(f"logit slope {slope.slope:.2f}, 95% CI {slope.slope_ci[0]:.2f}-{slope.slope_ci[1]:.2f}, " + f"p(slope=1) = {slope.slope_p_value:.1e}") +``` + +## Scores + +`brier_score` is the mean squared error `(p - y)^2`. A constant forecast of +0.5 scores 0.25. `log_loss` is the mean negative log-likelihood. It punishes +confident errors far more than the Brier score does. Forecasts are clipped to +`[eps, 1 - eps]` (default `1e-15`), because one market priced at exactly 0 or +1 that resolves the other way would otherwise make the loss infinite. Report +how many prices sit at 0 or 1 whenever you quote a log loss. + +Both accept `weights`, for example to weight markets by traded volume. + +## Brier decomposition + +Murphy (1973) splits the Brier score over forecast bins `k`, each with weight +share `W_k / W`, mean forecast `f_k` and observed frequency `o_k`: + +``` +brier = reliability - resolution + uncertainty +reliability = sum_k (W_k / W) (f_k - o_k)^2 lower is better +resolution = sum_k (W_k / W) (o_k - o)^2 higher is better +uncertainty = o (1 - o) set by the outcomes alone +``` + +This three-term identity is exact only when all forecasts in a bin are equal. +`brier_decomposition(bins=None)`, the default, uses one bin per distinct +forecast value, which suits prices quoted on a tick grid. With explicit `bins`, +forecasts vary within a bin, and two more terms keep the identity exact +(Stephenson, Coelho and Jolliffe, 2008): `within_bin_variance` and +`within_bin_covariance`, with +`brier = reliability - resolution + uncertainty + within_bin_variance - +within_bin_covariance`. `brier_skill_score` is `1 - brier / uncertainty`. + +## Reliability table + +`reliability_table` returns one row per bin: `n`, `n_clusters`, the summed +`weight`, `mean_forecast`, `observed_frequency`, `difference` (observed minus +forecast; negative means the bin is overpriced), and an interval for the +observed frequency. Bins are half-open `[lower, upper)`, and the last bin also +holds forecasts equal to 1. Empty bins remain in the table with null +statistics. + +The default edges, `PREDICTION_MARKET_BIN_EDGES`, are +`0, .01, .02, .05, .1, .2, ..., .9, .95, .98, .99, 1`. Market prices pile up +near 0 and 1, and the favorite-longshot bias is concentrated there. Equal-width +deciles would put every price below 10 cents into one bin. + +The interval method depends on the inputs and is recorded in `ci_method`: + +| Inputs | `ci_method` | Interval | +|---|---|---| +| no clusters, no weights | `wilson` | Wilson score interval | +| weights, no clusters | `wilson_effective_n` | Wilson on the Kish effective sample size `(sum w)^2 / sum w^2` | +| `clusters` given | `cluster_bootstrap` | Percentile interval from resampling whole clusters with replacement | + +Pass `clusters` whenever outcomes are dependent within a group. Several +markets on one election, or every strike of one price-range contract, resolve +together. A Wilson interval counts them as independent draws and can be too +narrow by a factor of up to `sqrt(markets per event)`. The bootstrap resamples +clusters across all bins at once. Replicates in which a bin is empty are left +out of that bin's interval, so read intervals for bins with few clusters with +care. `n_bootstrap` (default 1000) and `seed` (default 0, so repeated calls +match) control the resampling. + +## Calibration error summaries + +`expected_calibration_error` is `sum_k (W_k / W) |o_k - f_k|` over non-empty +bins, and `max_calibration_error` is the largest `|o_k - f_k|`. Both depend +on the bin edges, so report the edges with them. Sampling noise alone makes +the expected ECE positive for a perfectly calibrated forecast, and MCE is +often set by the thinnest bin. Check that bin's `n` and interval in the +reliability table. + +## Calibration slope and the favorite-longshot bias + +`calibration_slope` regresses outcomes on forecasts: + +- `spec="linear"`: `outcome = a + b * forecast`, weighted least squares (a + Mincer-Zarnowitz regression on probabilities). +- `spec="logit"`: `logit P(outcome) = a + b * logit(forecast)`, a logistic + regression (Cox's calibration slope). Forecasts are clipped to + `[eps, 1 - eps]` (default `eps=1e-3`) before the transform, because a price + of 0 or 1 has an infinite logit. + +A calibrated forecast has `b = 1` and `a = 0`. **`b > 1` is the +favorite-longshot bias**: outcomes rise faster than prices, so longshots win +less often than priced (overpriced) and favorites more often (underpriced). +`b < 1` is the reverse, the signature of an overconfident forecaster whose +probabilities are too extreme. + +Standard errors use the sandwich estimator with the small-sample factor +`G / (G - 1) * (n - 1) / (n - 2)`, where `G` is the number of clusters +(cluster-robust CR1) or `n` without clusters (HC1). Intervals use a t +distribution with `G - 1` degrees of freedom under clustering. `slope_p_value` +tests `b = 1`. + +A slope estimate is a summary. A bias confined to a few cents at either end +can leave `b` close to 1, so read the reliability table as well. And a +statistically significant `b` is not yet a trading edge: whether it survives +fees, spreads and capacity limits has to be checked separately. + +## References + +- Brier, G. W. (1950). "Verification of Forecasts Expressed in Terms of + Probability." *Monthly Weather Review*, 78(1), 1-3. +- Murphy, A. H. (1973). "A New Vector Partition of the Probability Score." + *Journal of Applied Meteorology*, 12(4), 595-600. +- Stephenson, D. B., Coelho, C. A. S., & Jolliffe, I. T. (2008). "Two Extra + Components in the Brier Score Decomposition." *Weather and Forecasting*, + 23(4), 752-757. +- Wilson, E. B. (1927). "Probable Inference, the Law of Succession, and + Statistical Inference." *Journal of the American Statistical Association*, + 22(158), 209-212. +- Cox, D. R. (1958). "Two Further Applications of a Model for Binary + Regression." *Biometrika*, 45(3/4), 562-565. +- Mincer, J., & Zarnowitz, V. (1969). "The Evaluation of Economic Forecasts." + In *Economic Forecasts and Expectations*, NBER, 3-46. +- Cameron, A. C., & Miller, D. L. (2015). "A Practitioner's Guide to + Cluster-Robust Inference." *Journal of Human Resources*, 50(2), 317-372. +- Snowberg, E., & Wolfers, J. (2010). "Explaining the Favorite-Long Shot Bias: + Is it Risk-Love or Misperceptions?" *Journal of Political Economy*, 118(4), + 723-746. diff --git a/docs/methods/index.md b/docs/methods/index.md index 2a36a1ed..9c10b9a4 100644 --- a/docs/methods/index.md +++ b/docs/methods/index.md @@ -28,6 +28,7 @@ The methods below address these problems with mathematical rigor. | [Deflated Sharpe Ratio](deflated-sharpe-ratio.md) | Selection bias from testing many strategies | `deflated_sharpe_ratio()` | Lopez de Prado et al. (2025) | | [CPCV](cpcv.md) | Backtest overfitting detection | `CombinatorialCV` | Lopez de Prado (2018) | | [HAC-adjusted IC](hac-ic.md) | Autocorrelation in IC significance testing | `compute_ic_hac_stats()` | Newey & West (1987) | +| [Probability Calibration](calibration.md) | Whether forecast probabilities or market prices match outcome frequencies | `reliability_table()`, `calibration_slope()` | Murphy (1973) | ## Methods by Category diff --git a/docs/reference/references.md b/docs/reference/references.md index 8a0e0518..4320d51b 100644 --- a/docs/reference/references.md +++ b/docs/reference/references.md @@ -577,6 +577,22 @@ High turnover means factor definitions change frequently, inducing trading costs "Testing and Valuing Dynamic Correlations for Asset Allocation." *Journal of Business & Economic Statistics*, 24(2), 238-253. +11. **Murphy, A. H. (1973)** + "A New Vector Partition of the Probability Score." + *Journal of Applied Meteorology*, 12(4), 595-600. + +12. **Stephenson, D. B., Coelho, C. A. S., & Jolliffe, I. T. (2008)** + "Two Extra Components in the Brier Score Decomposition." + *Weather and Forecasting*, 23(4), 752-757. + +13. **Cameron, A. C., & Miller, D. L. (2015)** + "A Practitioner's Guide to Cluster-Robust Inference." + *Journal of Human Resources*, 50(2), 317-372. + +14. **Snowberg, E., & Wolfers, J. (2010)** + "Explaining the Favorite-Long Shot Bias: Is it Risk-Love or Misperceptions?" + *Journal of Political Economy*, 118(4), 723-746. + ### Reference Implementations - DSR: https://github.com/zoonek/2025-sharpe-ratio diff --git a/mkdocs.yml b/mkdocs.yml index e4d5d4d1..f072e52e 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -126,6 +126,7 @@ nav: - Deflated Sharpe Ratio: methods/deflated-sharpe-ratio.md - Combinatorial Purged CV: methods/cpcv.md - HAC-adjusted IC: methods/hac-ic.md + - Probability Calibration: methods/calibration.md - Reference: - Architecture: reference/architecture.md - Academic References: reference/references.md diff --git a/src/ml4t/diagnostic/metrics/AGENTS.md b/src/ml4t/diagnostic/metrics/AGENTS.md index eaa466d2..4c2a37d8 100644 --- a/src/ml4t/diagnostic/metrics/AGENTS.md +++ b/src/ml4t/diagnostic/metrics/AGENTS.md @@ -12,6 +12,8 @@ Feature-level and signal-level metric functions used by diagnostics, selection, - `interactions.py` - H-statistic and SHAP interaction workflows - `feature_outcome.py` - combined feature-outcome analysis - `monotonicity.py` and `risk_adjusted.py` - ranking and portfolio-style metrics +- `calibration.py` - Brier score and decomposition, log loss, reliability table with + Wilson or cluster-bootstrap intervals, ECE/MCE, calibration slope ## Common Entry Points @@ -23,3 +25,4 @@ Feature-level and signal-level metric functions used by diagnostics, selection, - `compute_shap_importance` - `compute_h_statistic` - `analyze_feature_outcome` +- `reliability_table`, `brier_decomposition`, `calibration_slope` diff --git a/src/ml4t/diagnostic/metrics/__init__.py b/src/ml4t/diagnostic/metrics/__init__.py index fb8a2ccb..7e926161 100644 --- a/src/ml4t/diagnostic/metrics/__init__.py +++ b/src/ml4t/diagnostic/metrics/__init__.py @@ -7,6 +7,18 @@ from statsmodels.stats.sandwich_covariance import cov_hac from ml4t.diagnostic.metrics.basic import compute_forward_returns, hit_rate +from ml4t.diagnostic.metrics.calibration import ( + PREDICTION_MARKET_BIN_EDGES, + BrierDecomposition, + CalibrationSlope, + brier_decomposition, + brier_score, + calibration_slope, + expected_calibration_error, + log_loss, + max_calibration_error, + reliability_table, +) from ml4t.diagnostic.metrics.conditional import compute_conditional_ic from ml4t.diagnostic.metrics.feature_outcome import analyze_feature_outcome from ml4t.diagnostic.metrics.ic import ( @@ -87,4 +99,14 @@ "cross_sectional_auc_series", "compute_ic_uncertainty", "compute_auc_uncertainty", + "PREDICTION_MARKET_BIN_EDGES", + "BrierDecomposition", + "CalibrationSlope", + "brier_score", + "brier_decomposition", + "log_loss", + "reliability_table", + "expected_calibration_error", + "max_calibration_error", + "calibration_slope", ] diff --git a/src/ml4t/diagnostic/metrics/calibration.py b/src/ml4t/diagnostic/metrics/calibration.py new file mode 100644 index 00000000..d77939cf --- /dev/null +++ b/src/ml4t/diagnostic/metrics/calibration.py @@ -0,0 +1,807 @@ +"""Calibration diagnostics for probability forecasts of binary outcomes. + +A forecast is calibrated when events forecast with probability ``p`` occur with +frequency ``p``. The functions here score and decompose that property for any +probability forecast: model outputs, analyst probabilities, or prediction-market +prices read as implied probabilities. + +- ``brier_score`` and ``log_loss`` are proper scoring rules. +- ``brier_decomposition`` splits the Brier score into reliability, resolution + and uncertainty (Murphy, 1973), plus the two within-bin terms that make the + identity exact when forecasts vary inside a bin (Stephenson et al., 2008). +- ``reliability_table`` reports observed frequency against mean forecast per + probability bin, with a Wilson interval for independent outcomes or a + cluster-bootstrap interval when outcomes share an event. +- ``expected_calibration_error`` and ``max_calibration_error`` summarize the + table in one number. +- ``calibration_slope`` regresses outcomes on forecasts; a slope above one is + the favorite-longshot bias. + +All functions accept NumPy arrays, Polars or pandas Series, or sequences. +Forecasts must lie in ``[0, 1]`` and outcomes in ``{0, 1}``; missing values are +rejected rather than dropped so that a filtering decision stays with the caller. +""" + +from __future__ import annotations + +from collections.abc import Iterator, Sequence +from dataclasses import dataclass +from typing import Any, Literal + +import numpy as np +import pandas as pd +import polars as pl +import statsmodels.api as sm +from numpy.typing import NDArray +from scipy import sparse, stats + +ArrayLike = NDArray[Any] | pl.Series | pd.Series | Sequence[float] + +PREDICTION_MARKET_BIN_EDGES: tuple[float, ...] = ( + 0.0, + 0.01, + 0.02, + 0.05, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 0.95, + 0.98, + 0.99, + 1.0, +) +"""Default reliability bins: one-cent bins at the extremes, deciles in the middle. + +Prediction-market prices cluster near 0 and 1, and that is where the +favorite-longshot bias concentrates, so equal-width deciles would hide it in +the first and last bin. +""" + +_BOOTSTRAP_CHUNK_CELLS = 4_000_000 + + +# ============================================================================ +# Result types +# ============================================================================ + + +@dataclass(frozen=True) +class BrierDecomposition: + """Murphy decomposition of the Brier score. + + The identity + ``brier = reliability - resolution + uncertainty + within_bin_variance - + within_bin_covariance`` holds exactly. The two within-bin terms are zero + when every bin holds a single forecast value, which is the default binning + of ``brier_decomposition``. + + Attributes + ---------- + brier + Mean squared difference between forecast and outcome. + reliability + Weighted mean squared gap between each bin's mean forecast and its + observed frequency. Lower is better; zero is perfectly calibrated. + resolution + Weighted variance of bin observed frequencies around the overall base + rate. Higher is better: the forecast separates events from non-events. + uncertainty + ``base_rate * (1 - base_rate)``, the Brier score of always forecasting + the base rate. It depends on the outcomes only. + within_bin_variance + Weighted mean squared deviation of forecasts from their bin mean. + within_bin_covariance + Twice the weighted mean within-bin covariance of forecasts and outcomes. + base_rate + Weighted mean outcome. + n_bins + Number of non-empty bins used. + n + Number of observations. + """ + + brier: float + reliability: float + resolution: float + uncertainty: float + within_bin_variance: float + within_bin_covariance: float + base_rate: float + n_bins: int + n: int + + @property + def brier_skill_score(self) -> float: + """``1 - brier / uncertainty``; positive beats the base-rate forecast.""" + if self.uncertainty == 0.0: + return float("nan") + return 1.0 - self.brier / self.uncertainty + + +@dataclass(frozen=True) +class CalibrationSlope: + """Regression of outcomes on forecasts. + + ``outcome = intercept + slope * forecast`` for ``spec="linear"`` (a + Mincer-Zarnowitz regression on probabilities), or + ``logit P(outcome) = intercept + slope * logit(forecast)`` for + ``spec="logit"`` (Cox's calibration slope). A calibrated forecast has + ``slope = 1`` and ``intercept = 0``. + + ``slope > 1``: outcomes rise faster than forecasts, so low-probability + events (longshots) happen less often than priced and high-probability + events (favorites) more often. Longshots are overpriced and favorites + underpriced: the favorite-longshot bias. ``slope < 1`` is the reverse, the + pattern of an overconfident forecaster. + + Attributes + ---------- + intercept, slope + Point estimates. + intercept_se, slope_se + Sandwich standard errors with the small-sample factor + ``G / (G - 1) * (n - 1) / (n - 2)``: cluster-robust (CR1) with ``G`` + clusters when clusters are given, heteroskedasticity-robust (HC1, + ``G = n``) otherwise. + intercept_ci, slope_ci + Wald intervals at ``confidence``, using a t distribution with + ``n_clusters - 1`` degrees of freedom under clustering and ``n - 2`` + otherwise. + slope_p_value + Two-sided p-value for the null ``slope = 1``. + spec + ``"linear"`` or ``"logit"``. + se_type + ``"cluster"`` or ``"HC1"``. + confidence + Confidence level of the intervals. + n + Number of observations. + n_clusters + Number of distinct clusters, or ``n`` without clusters. + """ + + intercept: float + slope: float + intercept_se: float + slope_se: float + intercept_ci: tuple[float, float] + slope_ci: tuple[float, float] + slope_p_value: float + spec: str + se_type: str + confidence: float + n: int + n_clusters: int + + +# ============================================================================ +# Input handling +# ============================================================================ + + +def _to_numpy(values: ArrayLike | None) -> NDArray[Any] | None: + if values is None: + return None + if isinstance(values, pl.Series | pd.Series): + return values.to_numpy() + return np.asarray(values) + + +def _prepare( + forecasts: ArrayLike, + outcomes: ArrayLike, + weights: ArrayLike | None = None, +) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: + f_raw = _to_numpy(forecasts) + y_raw = _to_numpy(outcomes) + assert f_raw is not None and y_raw is not None + if f_raw.ndim != 1 or y_raw.ndim != 1: + raise ValueError("forecasts and outcomes must be one-dimensional") + if len(f_raw) != len(y_raw): + raise ValueError( + f"forecasts and outcomes must have the same length, got {len(f_raw)} and {len(y_raw)}" + ) + if len(f_raw) == 0: + raise ValueError("forecasts and outcomes must not be empty") + + f = f_raw.astype(np.float64) + y = y_raw.astype(np.float64) + if not np.all(np.isfinite(f)): + raise ValueError("forecasts contain NaN or infinite values") + if not np.all(np.isfinite(y)): + raise ValueError("outcomes contain NaN or infinite values") + if f.min() < 0.0 or f.max() > 1.0: + raise ValueError("forecasts must lie in [0, 1]") + if not np.all((y == 0.0) | (y == 1.0)): + raise ValueError("outcomes must be 0 or 1") + + w_raw = _to_numpy(weights) + if w_raw is None: + w = np.ones_like(f) + else: + if w_raw.shape != f.shape: + raise ValueError("weights must have the same length as forecasts") + w = w_raw.astype(np.float64) + if not np.all(np.isfinite(w)) or w.min() < 0.0: + raise ValueError("weights must be finite and non-negative") + if w.sum() <= 0.0: + raise ValueError("weights must not all be zero") + return f, y, w + + +def _cluster_codes(clusters: ArrayLike, n: int) -> tuple[NDArray[np.int64], int]: + raw = _to_numpy(clusters) + assert raw is not None + if raw.shape != (n,): + raise ValueError("clusters must have the same length as forecasts") + if pd.isna(raw).any(): + raise ValueError("clusters contain missing values") + codes, uniques = pd.factorize(raw, sort=False) + return codes.astype(np.int64), len(uniques) + + +def _validate_edges(bins: Sequence[float]) -> NDArray[np.float64]: + # Rounding snaps computed edges such as np.linspace(0, 1, 11)[3] == + # 0.30000000000000004 to the decimal a price of 0.3 is stored as, so a + # price on an edge lands in the bin above it as documented. + edges = np.round(np.asarray(bins, dtype=np.float64), 12) + if edges.ndim != 1 or len(edges) < 2: + raise ValueError("bins must contain at least two edges") + if not np.all(np.diff(edges) > 0): + raise ValueError("bins must be strictly increasing") + if edges[0] != 0.0 or edges[-1] != 1.0: + raise ValueError("bins must start at 0 and end at 1") + return edges + + +def _bin_index(f: NDArray[np.float64], edges: NDArray[np.float64]) -> NDArray[np.int64]: + """Assign forecasts to half-open bins ``[lo, hi)``; the last bin includes 1.""" + idx = np.searchsorted(edges, f, side="right") - 1 + return np.clip(idx, 0, len(edges) - 2).astype(np.int64) + + +def _z(confidence: float) -> float: + if not 0.0 < confidence < 1.0: + raise ValueError("confidence must be in (0, 1)") + return float(stats.norm.ppf(0.5 + confidence / 2.0)) + + +# ============================================================================ +# Scores +# ============================================================================ + + +def brier_score( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + weights: ArrayLike | None = None, +) -> float: + """Mean squared error of probability forecasts against binary outcomes. + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + weights + Optional non-negative observation weights. + + Returns + ------- + float + Weighted mean of ``(forecast - outcome) ** 2``; 0 is perfect, 0.25 is + the score of a constant 0.5 forecast. + """ + f, y, w = _prepare(forecasts, outcomes, weights) + return float(np.average((f - y) ** 2, weights=w)) + + +def log_loss( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + weights: ArrayLike | None = None, + eps: float = 1e-15, +) -> float: + """Mean negative log-likelihood of binary outcomes under the forecasts. + + Forecasts are clipped to ``[eps, 1 - eps]`` first. Without clipping, a + single market priced at 0 or 1 that resolves the other way makes the loss + infinite. With the default ``eps`` such an observation contributes about + 34.5, so report how many forecasts were clipped when prices of exactly 0 + or 1 occur in the data. + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + weights + Optional non-negative observation weights. + eps + Clipping bound, in ``(0, 0.5)``. + + Returns + ------- + float + Weighted mean of ``-(y log p + (1 - y) log(1 - p))`` in nats. + """ + if not 0.0 < eps < 0.5: + raise ValueError("eps must be in (0, 0.5)") + f, y, w = _prepare(forecasts, outcomes, weights) + p = np.clip(f, eps, 1.0 - eps) + losses = -(y * np.log(p) + (1.0 - y) * np.log1p(-p)) + return float(np.average(losses, weights=w)) + + +def brier_decomposition( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + bins: Sequence[float] | None = None, + weights: ArrayLike | None = None, +) -> BrierDecomposition: + """Decompose the Brier score into reliability, resolution and uncertainty. + + With ``bins=None`` every distinct forecast value is its own bin, which + gives Murphy's (1973) exact decomposition with zero within-bin terms. That + is the natural choice for prices quoted on a tick grid (prediction markets + quote in cents, so there are at most 101 values). For continuous model + probabilities pass bin edges; the two within-bin terms then absorb the + forecast variation inside each bin, keeping the identity exact + (Stephenson, Coelho and Jolliffe, 2008). + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + bins + Strictly increasing bin edges from 0 to 1, or ``None`` for one bin per + distinct forecast value. + weights + Optional non-negative observation weights. + + Returns + ------- + BrierDecomposition + Components satisfying ``brier = reliability - resolution + + uncertainty + within_bin_variance - within_bin_covariance``. + """ + f, y, w = _prepare(forecasts, outcomes, weights) + if bins is None: + _, idx = np.unique(f, return_inverse=True) + idx = idx.astype(np.int64) + else: + idx = _bin_index(f, _validate_edges(bins)) + + k = int(idx.max()) + 1 + total = w.sum() + w_k = np.bincount(idx, weights=w, minlength=k) + nonempty = w_k > 0 + f_bar = np.zeros(k) + o_bar = np.zeros(k) + f_bar[nonempty] = np.bincount(idx, weights=w * f, minlength=k)[nonempty] / w_k[nonempty] + o_bar[nonempty] = np.bincount(idx, weights=w * y, minlength=k)[nonempty] / w_k[nonempty] + base_rate = float((w * y).sum() / total) + + df = f - f_bar[idx] + dy = y - o_bar[idx] + return BrierDecomposition( + brier=float((w * (f - y) ** 2).sum() / total), + reliability=float((w_k * (f_bar - o_bar) ** 2).sum() / total), + resolution=float((w_k * (o_bar - base_rate) ** 2).sum() / total), + uncertainty=base_rate * (1.0 - base_rate), + within_bin_variance=float((w * df**2).sum() / total), + within_bin_covariance=float(2.0 * (w * df * dy).sum() / total), + base_rate=base_rate, + n_bins=int(nonempty.sum()), + n=len(f), + ) + + +# ============================================================================ +# Reliability +# ============================================================================ + + +def _wilson( + p_hat: NDArray[np.float64], n: NDArray[np.float64], z: float +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + """Vectorized Wilson score interval; ``n`` may be an effective sample size.""" + with np.errstate(invalid="ignore", divide="ignore"): + denom = 1.0 + z**2 / n + center = (p_hat + z**2 / (2.0 * n)) / denom + margin = z * np.sqrt((p_hat * (1.0 - p_hat) + z**2 / (4.0 * n)) / n) / denom + return center - margin, center + margin + + +def _cluster_bootstrap_counts( + n_clusters: int, + n_bootstrap: int, + rng: np.random.Generator, +) -> Iterator[NDArray[np.float64]]: + """Yield blocks of cluster multiplicities, one row per bootstrap replicate. + + Each row counts how often each cluster is drawn when ``n_clusters`` clusters + are sampled with replacement, so a statistic built from per-cluster sums is + recomputed for a replicate by a matrix product with the row. + """ + chunk = max(1, _BOOTSTRAP_CHUNK_CELLS // n_clusters) + done = 0 + while done < n_bootstrap: + b = min(chunk, n_bootstrap - done) + draws = rng.integers(0, n_clusters, size=(b, n_clusters)) + offsets = (np.arange(b) * n_clusters)[:, None] + counts = np.bincount((draws + offsets).ravel(), minlength=b * n_clusters) + yield counts.reshape(b, n_clusters).astype(np.float64) + done += b + + +def reliability_table( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + bins: Sequence[float] = PREDICTION_MARKET_BIN_EDGES, + clusters: ArrayLike | None = None, + weights: ArrayLike | None = None, + confidence: float = 0.95, + n_bootstrap: int = 1000, + seed: int | np.random.Generator | None = 0, +) -> pl.DataFrame: + """Observed outcome frequency against mean forecast, per probability bin. + + Bins are half-open ``[lower, upper)``; the last bin also holds forecasts + equal to 1. Every bin appears in the output, and statistics of an empty bin + are null. + + The interval for the observed frequency depends on the inputs: + + - no clusters, no weights: Wilson score interval (``ci_method="wilson"``); + - weights, no clusters: Wilson interval on the Kish effective sample size + ``(sum w)^2 / sum w^2`` (``"wilson_effective_n"``); + - clusters: percentile interval from a cluster bootstrap that resamples + whole clusters with replacement across all bins (``"cluster_bootstrap"``). + Use this when outcomes are correlated within a group, for example + several markets on one event. A Wilson interval treats them as + independent and is too narrow. Bootstrap replicates in which a bin is + empty are left out of that bin's interval. + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + bins + Strictly increasing bin edges from 0 to 1. The default, + ``PREDICTION_MARKET_BIN_EDGES``, is finer near 0 and 1. + clusters + Optional cluster label per observation (any hashable dtype). + weights + Optional non-negative observation weights. + confidence + Confidence level of the interval. + n_bootstrap + Number of cluster-bootstrap replicates; ignored without clusters. + seed + Seed or Generator for the cluster bootstrap. The default makes repeated + calls on the same data return the same interval; ``None`` draws a fresh + stream. + + Returns + ------- + pl.DataFrame + One row per bin with columns ``bin_lower``, ``bin_upper``, ``n``, + ``n_clusters``, ``weight``, ``mean_forecast``, ``observed_frequency``, + ``difference`` (observed minus forecast; negative means the bin is + overpriced), ``ci_lower``, ``ci_upper`` and ``ci_method``. + """ + f, y, w = _prepare(forecasts, outcomes, weights) + edges = _validate_edges(bins) + z = _z(confidence) + idx = _bin_index(f, edges) + k = len(edges) - 1 + + n_k = np.bincount(idx, minlength=k) + w_k = np.bincount(idx, weights=w, minlength=k) + wy_k = np.bincount(idx, weights=w * y, minlength=k) + wf_k = np.bincount(idx, weights=w * f, minlength=k) + nonempty = w_k > 0 + with np.errstate(invalid="ignore", divide="ignore"): + mean_f = np.where(nonempty, wf_k / w_k, np.nan) + obs = np.where(nonempty, wy_k / w_k, np.nan) + + if clusters is None: + nclu_k = n_k.copy() + if weights is None: + n_eff = n_k.astype(np.float64) + method = "wilson" + else: + w2_k = np.bincount(idx, weights=w**2, minlength=k) + with np.errstate(invalid="ignore", divide="ignore"): + n_eff = np.where(w2_k > 0, w_k**2 / w2_k, 0.0) + method = "wilson_effective_n" + lo, hi = _wilson(obs, n_eff, z) + else: + if isinstance(n_bootstrap, bool) or not isinstance(n_bootstrap, int) or n_bootstrap < 2: + raise ValueError("n_bootstrap must be an integer of at least 2") + codes, g = _cluster_codes(clusters, len(f)) + nclu_k = np.bincount(np.unique(codes * k + idx) % k, minlength=k) + rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed) + # Per-cluster, per-bin sums; a replicate's bin sums are counts @ sums. + s_w = sparse.csr_matrix((w, (codes, idx)), shape=(g, k)) + s_wy = sparse.csr_matrix((w * y, (codes, idx)), shape=(g, k)) + replicates = [] + for counts in _cluster_bootstrap_counts(g, n_bootstrap, rng): + bw = np.asarray((s_w.T @ counts.T).T) + bwy = np.asarray((s_wy.T @ counts.T).T) + with np.errstate(invalid="ignore", divide="ignore"): + replicates.append(np.where(bw > 0, bwy / bw, np.nan)) + boot = np.vstack(replicates) + alpha = 1.0 - confidence + lo = np.full(k, np.nan) + hi = np.full(k, np.nan) + for j in np.flatnonzero(nonempty): + col = boot[:, j] + col = col[~np.isnan(col)] + if col.size: + lo[j], hi[j] = np.quantile(col, [alpha / 2.0, 1.0 - alpha / 2.0]) + method = "cluster_bootstrap" + + def _nullable(values: NDArray[np.float64]) -> list[float | None]: + return [float(v) if m else None for v, m in zip(values, nonempty, strict=True)] + + return pl.DataFrame( + { + "bin_lower": edges[:-1], + "bin_upper": edges[1:], + "n": n_k.astype(np.int64), + "n_clusters": nclu_k.astype(np.int64), + "weight": w_k, + "mean_forecast": _nullable(mean_f), + "observed_frequency": _nullable(obs), + "difference": _nullable(obs - mean_f), + "ci_lower": _nullable(lo), + "ci_upper": _nullable(hi), + "ci_method": [method] * k, + }, + schema={ + "bin_lower": pl.Float64, + "bin_upper": pl.Float64, + "n": pl.Int64, + "n_clusters": pl.Int64, + "weight": pl.Float64, + "mean_forecast": pl.Float64, + "observed_frequency": pl.Float64, + "difference": pl.Float64, + "ci_lower": pl.Float64, + "ci_upper": pl.Float64, + "ci_method": pl.String, + }, + ) + + +def _bin_gaps( + forecasts: ArrayLike, + outcomes: ArrayLike, + bins: Sequence[float], + weights: ArrayLike | None, +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + f, y, w = _prepare(forecasts, outcomes, weights) + edges = _validate_edges(bins) + idx = _bin_index(f, edges) + k = len(edges) - 1 + w_k = np.bincount(idx, weights=w, minlength=k) + nonempty = w_k > 0 + w_k = w_k[nonempty] + f_bar = np.bincount(idx, weights=w * f, minlength=k)[nonempty] / w_k + o_bar = np.bincount(idx, weights=w * y, minlength=k)[nonempty] / w_k + return np.abs(o_bar - f_bar), w_k / w_k.sum() + + +def expected_calibration_error( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + bins: Sequence[float] = PREDICTION_MARKET_BIN_EDGES, + weights: ArrayLike | None = None, +) -> float: + """Weight-averaged absolute gap between observed frequency and mean forecast. + + ``ECE = sum_k (W_k / W) * |observed_k - forecast_k|`` over non-empty bins, + using the same bins as ``reliability_table``. The value depends on the bin + edges, so report them with it. Sampling noise inflates ECE in thin bins: + a perfectly calibrated forecast has positive expected ECE. + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + bins + Strictly increasing bin edges from 0 to 1. + weights + Optional non-negative observation weights. + + Returns + ------- + float + Expected calibration error in probability units. + """ + gaps, shares = _bin_gaps(forecasts, outcomes, bins, weights) + return float((gaps * shares).sum()) + + +def max_calibration_error( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + bins: Sequence[float] = PREDICTION_MARKET_BIN_EDGES, + weights: ArrayLike | None = None, +) -> float: + """Largest absolute gap between observed frequency and mean forecast. + + Taken over non-empty bins. A bin with a handful of observations can + dominate this number; check ``n`` and the interval in + ``reliability_table`` for the bin that attains it. + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + bins + Strictly increasing bin edges from 0 to 1. + weights + Optional non-negative observation weights. + + Returns + ------- + float + Maximum calibration error in probability units. + """ + gaps, _ = _bin_gaps(forecasts, outcomes, bins, weights) + return float(gaps.max()) + + +# ============================================================================ +# Calibration slope +# ============================================================================ + + +def calibration_slope( + forecasts: ArrayLike, + outcomes: ArrayLike, + *, + spec: Literal["linear", "logit"] = "linear", + clusters: ArrayLike | None = None, + weights: ArrayLike | None = None, + confidence: float = 0.95, + eps: float = 1e-3, +) -> CalibrationSlope: + """Regress outcomes on forecasts to measure the favorite-longshot bias. + + ``spec="linear"`` fits ``outcome = a + b * forecast`` by (weighted) least + squares, a Mincer-Zarnowitz regression on probabilities. ``spec="logit"`` + fits the logistic regression ``logit P(outcome) = a + b * logit(forecast)``, + Cox's (1958) calibration slope, which respects the ``[0, 1]`` range and + weights the tails more than the linear form does. + + Interpretation, for either specification: + + - ``b = 1`` and ``a = 0``: calibrated. + - ``b > 1``: favorite-longshot bias. Low-priced outcomes (longshots) win + less often than their price implies and high-priced outcomes + (favorites) more often, so longshots are overpriced and favorites + underpriced. + - ``b < 1``: reverse bias; forecasts are too extreme (overconfident). + + Standard errors are cluster-robust (CR1) when ``clusters`` is given and + HC1 otherwise. Pass clusters whenever outcomes are correlated within a + group, such as markets that resolve on one event; ignoring that + correlation understates the standard errors. + + Parameters + ---------- + forecasts + Probabilities in ``[0, 1]``. + outcomes + Realized outcomes, 0 or 1. + spec + ``"linear"`` or ``"logit"``. + clusters + Optional cluster label per observation. + weights + Optional non-negative observation weights. + confidence + Confidence level of the Wald intervals. + eps + For ``spec="logit"``, forecasts are clipped to ``[eps, 1 - eps]`` + before the logit transform. A price of 0 or 1 has an infinite logit, + and very small ``eps`` gives such observations extreme leverage. + Ignored for ``spec="linear"``. + + Returns + ------- + CalibrationSlope + Coefficients, standard errors, intervals and the p-value for + ``slope = 1``. + """ + if spec not in ("linear", "logit"): + raise ValueError("spec must be 'linear' or 'logit'") + if not 0.0 < eps < 0.5: + raise ValueError("eps must be in (0, 0.5)") + _z(confidence) + f, y, w = _prepare(forecasts, outcomes, weights) + + if spec == "logit": + p = np.clip(f, eps, 1.0 - eps) + x = np.log(p) - np.log1p(-p) + else: + x = f + if np.ptp(x) == 0.0: + raise ValueError("forecasts are constant; the slope is not identified") + exog = np.column_stack([np.ones_like(x), x]) + + n = len(f) + if spec == "linear": + sw = np.sqrt(w) + params = np.linalg.lstsq(exog * sw[:, None], y * sw, rcond=None)[0] + mu = exog @ params + hessian = exog.T @ (exog * w[:, None]) + else: + glm = sm.GLM(y, exog, family=sm.families.Binomial(), var_weights=w).fit() + params = np.asarray(glm.params, dtype=np.float64) + mu = 1.0 / (1.0 + np.exp(-(exog @ params))) + hessian = exog.T @ (exog * (w * mu * (1.0 - mu))[:, None]) + scores = exog * (w * (y - mu))[:, None] + + if clusters is None: + n_groups = n + meat_scores = scores + se_type = "HC1" + else: + codes, n_groups = _cluster_codes(clusters, n) + if n_groups < 2: + raise ValueError("cluster-robust standard errors need at least two clusters") + meat_scores = np.column_stack( + [np.bincount(codes, weights=scores[:, j], minlength=n_groups) for j in range(2)] + ) + se_type = "cluster" + bread = np.linalg.inv(hessian) + correction = n_groups / (n_groups - 1) * (n - 1) / (n - 2) + cov = correction * bread @ (meat_scores.T @ meat_scores) @ bread + bse = np.sqrt(np.diag(cov)) + + t_dist = stats.t(df=n_groups - 1 if clusters is not None else n - 2) + q = float(t_dist.ppf(0.5 + confidence / 2.0)) + ci = np.column_stack([params - q * bse, params + q * bse]) + slope_p_value = float(2.0 * t_dist.sf(abs((params[1] - 1.0) / bse[1]))) + return CalibrationSlope( + intercept=float(params[0]), + slope=float(params[1]), + intercept_se=float(bse[0]), + slope_se=float(bse[1]), + intercept_ci=(float(ci[0, 0]), float(ci[0, 1])), + slope_ci=(float(ci[1, 0]), float(ci[1, 1])), + slope_p_value=slope_p_value, + spec=spec, + se_type=se_type, + confidence=confidence, + n=n, + n_clusters=int(n_groups), + ) diff --git a/tests/metrics/test_calibration.py b/tests/metrics/test_calibration.py new file mode 100644 index 00000000..669bed4c --- /dev/null +++ b/tests/metrics/test_calibration.py @@ -0,0 +1,421 @@ +"""Tests for calibration diagnostics of probability forecasts.""" + +from __future__ import annotations + +import math + +import numpy as np +import pandas as pd +import polars as pl +import pytest + +from ml4t.diagnostic.evaluation.binary_metrics import wilson_score_interval +from ml4t.diagnostic.metrics import ( + PREDICTION_MARKET_BIN_EDGES, + brier_decomposition, + brier_score, + calibration_slope, + expected_calibration_error, + log_loss, + max_calibration_error, + reliability_table, +) + +DECILES = tuple(np.linspace(0.0, 1.0, 11)) + + +def _logit(p: np.ndarray) -> np.ndarray: + return np.log(p) - np.log1p(-p) + + +def _expit(x: np.ndarray) -> np.ndarray: + return 1.0 / (1.0 + np.exp(-x)) + + +def _market_prices(n: int, rng: np.random.Generator) -> np.ndarray: + """Cent-quoted prices spread across [0.01, 0.99], denser near the extremes.""" + raw = rng.beta(0.6, 0.6, size=n) + return np.clip(np.round(raw, 2), 0.01, 0.99) + + +def _favorite_longshot(p: np.ndarray, b: float) -> np.ndarray: + """True win probability when prices are compressed toward 0.5 by factor b.""" + return _expit(b * _logit(p)) + + +# ============================================================================ +# Scores and Brier decomposition +# ============================================================================ + + +def test_brier_and_log_loss_match_hand_computed_values() -> None: + """Guards against a wrong loss formula or a sign error in the log terms.""" + f = [0.2, 0.2, 0.8, 0.8] + y = [0, 1, 1, 1] + assert brier_score(f, y) == pytest.approx((0.04 + 0.64 + 0.04 + 0.04) / 4) + expected_ll = -(math.log(0.8) + math.log(0.2) + 2 * math.log(0.8)) / 4 + assert log_loss(f, y) == pytest.approx(expected_ll) + + +def test_log_loss_clips_certain_wrong_forecasts_to_a_finite_value() -> None: + """A market priced at 0 that resolves YES must not make the loss infinite.""" + value = log_loss([0.0, 1.0], [1, 1]) + assert math.isfinite(value) + assert value == pytest.approx(-math.log(1e-15) / 2, rel=1e-6) + expected = -(math.log(0.01) + math.log(0.99)) / 2 + assert log_loss([0.0, 1.0], [1, 1], eps=0.01) == pytest.approx(expected) + + +def test_brier_decomposition_matches_hand_computed_components() -> None: + """Guards against swapping reliability and resolution or a wrong base rate.""" + d = brier_decomposition([0.2, 0.2, 0.8, 0.8], [0, 1, 1, 1]) + assert d.reliability == pytest.approx((2 * 0.3**2 + 2 * 0.2**2) / 4) + assert d.resolution == pytest.approx((2 * 0.25**2 + 2 * 0.25**2) / 4) + assert d.uncertainty == pytest.approx(0.75 * 0.25) + assert d.brier == pytest.approx(0.19) + assert d.n_bins == 2 + assert d.within_bin_variance == 0.0 + assert d.within_bin_covariance == 0.0 + + +@pytest.mark.parametrize("bins", [None, DECILES, PREDICTION_MARKET_BIN_EDGES]) +def test_brier_decomposition_identity_is_exact(bins: tuple[float, ...] | None) -> None: + """The components must add back to the Brier score for any binning. + + With continuous forecasts inside wide bins the classic three-term identity + does not hold; dropping either within-bin term breaks this test. + """ + rng = np.random.default_rng(7) + f = rng.uniform(size=5_000) + y = (rng.uniform(size=5_000) < _favorite_longshot(f, 1.4)).astype(int) + w = rng.uniform(0.5, 2.0, size=5_000) + for weights in (None, w): + d = brier_decomposition(f, y, bins=bins, weights=weights) + recomposed = ( + d.reliability + - d.resolution + + d.uncertainty + + d.within_bin_variance + - d.within_bin_covariance + ) + assert recomposed == pytest.approx(d.brier, abs=1e-12) + assert d.brier == pytest.approx(brier_score(f, y, weights=weights)) + if bins is not None: + assert d.within_bin_variance > 0.0 + + +def test_brier_reliability_separates_calibrated_from_distorted_prices() -> None: + """Reliability must be near zero for calibrated prices and not for distorted ones.""" + rng = np.random.default_rng(11) + p = _market_prices(200_000, rng) + u = rng.uniform(size=p.size) + calibrated = brier_decomposition(p, (u < p).astype(int)) + distorted = brier_decomposition(p, (u < _favorite_longshot(p, 1.5)).astype(int)) + assert calibrated.reliability < 2e-4 + assert distorted.reliability > 10 * calibrated.reliability + + +# ============================================================================ +# Reliability table +# ============================================================================ + + +def test_reliability_table_recovers_planted_favorite_longshot_bias() -> None: + """Each bin's interval must cover the planted true frequency. + + Longshot bins must show observed below price and favorite bins above it, + with the price itself outside the interval in the tails. + """ + rng = np.random.default_rng(3) + p = _market_prices(300_000, rng) + q = _favorite_longshot(p, 1.5) + y = (rng.uniform(size=p.size) < q).astype(int) + + table = reliability_table(p, y) + edges = np.asarray(PREDICTION_MARKET_BIN_EDGES) + idx = np.clip(np.searchsorted(edges, p, side="right") - 1, 0, len(edges) - 2) + true_freq = np.bincount(idx, weights=q, minlength=len(edges) - 1) / np.bincount( + idx, minlength=len(edges) - 1 + ).clip(1) + + rows = table.filter(pl.col("n") > 0) + assert rows["ci_method"].unique().to_list() == ["wilson"] + covered = 0 + for row in rows.iter_rows(named=True): + j = int(np.searchsorted(edges, row["bin_lower"])) + covered += row["ci_lower"] <= true_freq[j] <= row["ci_upper"] + assert covered >= len(rows) - 1 + + low = rows.filter(pl.col("bin_upper") <= 0.1) + high = rows.filter(pl.col("bin_lower") >= 0.9) + assert (low["difference"] < 0).all() + assert (high["difference"] > 0).all() + assert (low["ci_upper"] < low["mean_forecast"]).all() + assert (high["ci_lower"] > high["mean_forecast"]).all() + + +def test_reliability_table_bins_are_half_open_and_include_one() -> None: + """A price on an edge belongs to the bin above it; a price of 1 to the last bin.""" + table = reliability_table([0.0, 0.01, 0.019, 0.02, 0.99, 1.0], [0, 0, 0, 1, 1, 1]) + counts = dict(zip(table["bin_lower"].to_list(), table["n"].to_list(), strict=True)) + assert counts[0.0] == 1 + assert counts[0.01] == 2 + assert counts[0.02] == 1 + assert counts[0.99] == 2 + assert table["n"].sum() == 6 + + +def test_computed_bin_edges_assign_edge_prices_to_the_upper_bin() -> None: + """np.linspace gives 0.30000000000000004; a price of 0.3 must still go to [0.3, 0.4).""" + assert DECILES[3] != 0.3 + table = reliability_table([0.3, 0.7], [0, 1], bins=DECILES) + occupied = table.filter(pl.col("n") > 0)["bin_lower"].to_list() + assert occupied == pytest.approx([0.3, 0.7]) + + +def test_reliability_table_empty_bins_are_null() -> None: + """Empty bins stay in the table with null statistics, not NaN or zeros.""" + table = reliability_table([0.3, 0.35], [0, 1], bins=DECILES) + assert table.height == 10 + empty = table.filter(pl.col("n") == 0) + assert empty.height == 9 + assert empty["observed_frequency"].null_count() == 9 + assert empty["ci_lower"].null_count() == 9 + + +def test_reliability_wilson_interval_matches_library_wilson() -> None: + """The vectorized interval must agree with ``wilson_score_interval``.""" + rng = np.random.default_rng(5) + p = rng.uniform(size=2_000) + y = (rng.uniform(size=2_000) < p).astype(int) + table = reliability_table(p, y, bins=DECILES, confidence=0.9) + for row in table.iter_rows(named=True): + k = round(row["observed_frequency"] * row["n"]) + lo, hi = wilson_score_interval(k, row["n"], confidence=0.9) + assert row["ci_lower"] == pytest.approx(lo) + assert row["ci_upper"] == pytest.approx(hi) + + +def test_cluster_bootstrap_widens_when_outcomes_are_shared_within_clusters() -> None: + """Perfectly correlated clusters of size m carry ~1/m of the information. + + Resampling observations instead of clusters, or ignoring the cluster + ids, gives intervals as narrow as the independent case and fails here. + """ + rng = np.random.default_rng(21) + n_clusters, size = 600, 9 + cluster_p = rng.uniform(0.05, 0.95, size=n_clusters) + p = np.repeat(cluster_p, size) + clusters = np.repeat(np.arange(n_clusters), size) + shared = np.repeat((rng.uniform(size=n_clusters) < cluster_p).astype(int), size) + independent = (rng.uniform(size=p.size) < p).astype(int) + + def widths(y: np.ndarray) -> np.ndarray: + t = reliability_table(p, y, bins=DECILES, clusters=clusters, n_bootstrap=500) + assert t["ci_method"].unique().to_list() == ["cluster_bootstrap"] + return (t["ci_upper"] - t["ci_lower"]).to_numpy() + + ratio = np.median(widths(shared) / widths(independent)) + assert ratio > 0.7 * math.sqrt(size) + + +def test_cluster_bootstrap_with_singleton_clusters_matches_wilson_width() -> None: + """With one observation per cluster the bootstrap must agree with Wilson.""" + rng = np.random.default_rng(8) + p = rng.uniform(0.1, 0.9, size=20_000) + y = (rng.uniform(size=p.size) < p).astype(int) + wilson = reliability_table(p, y, bins=DECILES) + boot = reliability_table(p, y, bins=DECILES, clusters=np.arange(p.size), n_bootstrap=400) + ratio = (boot["ci_upper"] - boot["ci_lower"]) / (wilson["ci_upper"] - wilson["ci_lower"]) + assert ratio.drop_nulls().is_between(0.8, 1.2).all() + assert boot["n_clusters"].to_list() == boot["n"].to_list() + + +def test_cluster_bootstrap_counts_clusters_per_bin_and_is_seeded() -> None: + """``n_clusters`` counts distinct clusters per bin; the seed fixes the interval.""" + p = np.array([0.15, 0.15, 0.15, 0.55, 0.55, 0.85]) + y = np.array([0, 1, 0, 1, 1, 1]) + clusters = pl.Series(["a", "a", "b", "c", "c", "a"]) + t1 = reliability_table(p, y, bins=DECILES, clusters=clusters, n_bootstrap=200, seed=4) + t2 = reliability_table(p, y, bins=DECILES, clusters=clusters, n_bootstrap=200, seed=4) + by_bin = dict(zip(t1["bin_lower"].round(2).to_list(), t1["n_clusters"].to_list(), strict=True)) + assert by_bin[0.1] == 2 + assert by_bin[0.5] == 1 + assert by_bin[0.8] == 1 + assert t1.equals(t2) + + rng = np.random.default_rng(0) + big_p = rng.uniform(size=3_000) + big_y = (rng.uniform(size=3_000) < big_p).astype(int) + big_c = rng.integers(0, 300, size=3_000) + a = reliability_table(big_p, big_y, bins=DECILES, clusters=big_c, n_bootstrap=200, seed=1) + b = reliability_table(big_p, big_y, bins=DECILES, clusters=big_c, n_bootstrap=200, seed=2) + assert not a["ci_lower"].equals(b["ci_lower"]) + + +def test_weights_equal_duplicated_observations() -> None: + """An observation with weight 2 must count exactly like two copies of it.""" + rng = np.random.default_rng(9) + p = rng.uniform(size=500) + y = (rng.uniform(size=500) < p).astype(int) + w = rng.integers(1, 4, size=500) + p_dup, y_dup = np.repeat(p, w), np.repeat(y, w) + + assert brier_score(p, y, weights=w) == pytest.approx(brier_score(p_dup, y_dup)) + assert log_loss(p, y, weights=w) == pytest.approx(log_loss(p_dup, y_dup)) + weighted = reliability_table(p, y, bins=DECILES, weights=w) + duplicated = reliability_table(p_dup, y_dup, bins=DECILES) + for column in ("mean_forecast", "observed_frequency", "difference"): + assert np.allclose(weighted[column].to_numpy(), duplicated[column].to_numpy()) + assert weighted["ci_method"][0] == "wilson_effective_n" + for spec in ("linear", "logit"): + a = calibration_slope(p, y, spec=spec, weights=w) + b = calibration_slope(p_dup, y_dup, spec=spec) + assert a.slope == pytest.approx(b.slope, rel=1e-6) + assert a.intercept == pytest.approx(b.intercept, abs=1e-6) + + +# ============================================================================ +# Calibration error summaries +# ============================================================================ + + +def test_ece_and_mce_match_hand_computed_values() -> None: + """Guards against an unweighted average over bins or including empty bins.""" + f = [0.15, 0.15, 0.15, 0.15, 0.85] + y = [0, 0, 0, 1, 0] + assert expected_calibration_error(f, y, bins=DECILES) == pytest.approx(0.8 * 0.1 + 0.2 * 0.85) + assert max_calibration_error(f, y, bins=DECILES) == pytest.approx(0.85) + + +def test_ece_is_small_when_calibrated_and_large_when_distorted() -> None: + """ECE must detect a planted distortion that calibrated data does not show.""" + rng = np.random.default_rng(13) + p = _market_prices(200_000, rng) + u = rng.uniform(size=p.size) + assert expected_calibration_error(p, (u < p).astype(int)) < 0.005 + assert expected_calibration_error(p, (u < _favorite_longshot(p, 1.5)).astype(int)) > 0.02 + + +# ============================================================================ +# Calibration slope +# ============================================================================ + + +@pytest.mark.parametrize("spec", ["linear", "logit"]) +def test_calibration_slope_is_one_for_calibrated_prices(spec: str) -> None: + """Calibrated prices must give a slope whose interval covers 1.""" + rng = np.random.default_rng(17) + p = _market_prices(100_000, rng) + y = (rng.uniform(size=p.size) < p).astype(int) + result = calibration_slope(p, y, spec=spec) + assert result.slope_ci[0] <= 1.0 <= result.slope_ci[1] + assert result.intercept_ci[0] <= 0.0 <= result.intercept_ci[1] + assert result.slope_p_value > 0.05 + assert result.slope == pytest.approx(1.0, abs=0.03) + + +def test_logit_calibration_slope_recovers_planted_slope() -> None: + """The logit slope must recover b in logit(q) = b * logit(p) and reject b = 1.""" + rng = np.random.default_rng(19) + p = _market_prices(100_000, rng) + y = (rng.uniform(size=p.size) < _favorite_longshot(p, 1.3)).astype(int) + result = calibration_slope(p, y, spec="logit") + assert result.slope_ci[0] <= 1.3 <= result.slope_ci[1] + assert result.slope_ci[0] > 1.0 + assert result.slope_p_value < 1e-6 + + +def test_linear_calibration_slope_signals_favorite_longshot_bias() -> None: + """Overpriced longshots and underpriced favorites must give b > 1, a < 0.""" + rng = np.random.default_rng(23) + p = _market_prices(100_000, rng) + y = (rng.uniform(size=p.size) < _favorite_longshot(p, 1.5)).astype(int) + result = calibration_slope(p, y, spec="linear") + assert result.slope_ci[0] > 1.0 + assert result.intercept_ci[1] < 0.0 + + +def test_calibration_slope_matches_statsmodels_robust_covariance() -> None: + """The hand-written sandwich must equal statsmodels' HC1 and CR1 estimates.""" + sm = pytest.importorskip("statsmodels.api") + rng = np.random.default_rng(29) + p = rng.uniform(0.02, 0.98, size=4_000) + y = (rng.uniform(size=p.size) < p).astype(int) + groups = rng.integers(0, 200, size=p.size) + exog = sm.add_constant(p) + + hc1 = sm.OLS(y, exog).fit(cov_type="HC1") + ours = calibration_slope(p, y) + assert [ours.intercept_se, ours.slope_se] == pytest.approx(list(hc1.bse), rel=1e-8) + + cr1 = sm.OLS(y, exog).fit(cov_type="cluster", cov_kwds={"groups": groups}) + ours = calibration_slope(p, y, clusters=groups) + assert [ours.intercept_se, ours.slope_se] == pytest.approx(list(cr1.bse), rel=1e-8) + assert ours.n_clusters == 200 + + logit_exog = sm.add_constant(np.log(p / (1 - p))) + glm = sm.GLM(y, logit_exog, family=sm.families.Binomial()).fit( + cov_type="cluster", cov_kwds={"groups": groups} + ) + ours = calibration_slope(p, y, spec="logit", clusters=groups) + assert [ours.intercept, ours.slope] == pytest.approx(list(glm.params), rel=1e-6) + assert [ours.intercept_se, ours.slope_se] == pytest.approx(list(glm.bse), rel=1e-6) + + +def test_cluster_robust_slope_se_widens_with_shared_outcomes() -> None: + """Ignoring within-event correlation understates the slope standard error.""" + rng = np.random.default_rng(31) + n_clusters, size = 2_000, 9 + cluster_p = rng.uniform(0.05, 0.95, size=n_clusters) + p = np.repeat(cluster_p, size) + clusters = np.repeat(np.arange(n_clusters), size) + y = np.repeat((rng.uniform(size=n_clusters) < cluster_p).astype(int), size) + clustered = calibration_slope(p, y, clusters=clusters) + naive = calibration_slope(p, y) + assert clustered.se_type == "cluster" + assert naive.se_type == "HC1" + assert clustered.slope_se / naive.slope_se > 0.8 * math.sqrt(size) + + +# ============================================================================ +# Inputs +# ============================================================================ + + +def test_accepts_polars_and_pandas_inputs() -> None: + """Series inputs must give the same answer as NumPy arrays.""" + f = np.array([0.1, 0.4, 0.6, 0.9]) + y = np.array([0, 0, 1, 1]) + expected = brier_score(f, y) + assert brier_score(pl.Series(f), pl.Series(y)) == pytest.approx(expected) + assert brier_score(pd.Series(f), pd.Series(y.astype(bool))) == pytest.approx(expected) + + +@pytest.mark.parametrize( + ("forecasts", "outcomes", "message"), + [ + ([0.1, np.nan], [0, 1], "NaN"), + ([0.1, 1.2], [0, 1], r"\[0, 1\]"), + ([0.1, 0.2], [0, 2], "0 or 1"), + ([0.1, 0.2], [0], "same length"), + ([], [], "empty"), + ], +) +def test_rejects_invalid_inputs(forecasts: list[float], outcomes: list[int], message: str) -> None: + """Invalid probabilities or outcomes must fail loudly instead of being dropped.""" + with pytest.raises(ValueError, match=message): + brier_score(forecasts, outcomes) + + +@pytest.mark.parametrize("bins", [(0.0, 0.5, 0.4, 1.0), (0.1, 1.0), (0.0, 0.9), (0.5,)]) +def test_rejects_invalid_bin_edges(bins: tuple[float, ...]) -> None: + """Bins must cover [0, 1] with strictly increasing edges.""" + with pytest.raises(ValueError, match="bins"): + reliability_table([0.2, 0.7], [0, 1], bins=bins) + + +def test_calibration_slope_rejects_constant_forecasts() -> None: + """A constant forecast does not identify a slope.""" + with pytest.raises(ValueError, match="constant"): + calibration_slope([0.5, 0.5, 0.5], [0, 1, 1]) diff --git a/tests/metrics/test_public_metrics_namespace.py b/tests/metrics/test_public_metrics_namespace.py index 4a2f4e52..bcf5c4a0 100644 --- a/tests/metrics/test_public_metrics_namespace.py +++ b/tests/metrics/test_public_metrics_namespace.py @@ -30,6 +30,13 @@ def test_metrics_package_exports_core_metric_groups(): "compute_fold_percentiles", "QuantileProfile", "quantile_profile", + "brier_score", + "brier_decomposition", + "log_loss", + "reliability_table", + "expected_calibration_error", + "max_calibration_error", + "calibration_slope", } assert expected.issubset(set(metrics.__all__)) diff --git a/uv.lock b/uv.lock index 4fe4642c..35654d7d 100644 --- a/uv.lock +++ b/uv.lock @@ -3026,11 +3026,11 @@ wheels = [ [[package]] name = "pypdf" -version = "6.16.2" +version = "6.19.0" source = { registry = "https://pypi.org/simple" } -sdist = { url = "https://files.pythonhosted.org/packages/44/66/54212e75406afd9f3e933d0dda23072f6aecc55c5a273077dc2e0b028b23/pypdf-6.16.2.tar.gz", hash = "sha256:595647f6191de6f402cfde1d0c455d6cbccbd509aac32b34783009c032de5d6e", size = 7008996, upload-time = "2026-08-23T13:50:07.135Z" } +sdist = { url = "https://files.pythonhosted.org/packages/1f/ac/63d71aaedb59acbcdef491e6ca6469165e3771c9c74358204818fd9bc5a6/pypdf-6.19.0.tar.gz", hash = "sha256:bbc43aca292369ccc6cbc8a921991ecf2538a3587ab5a116eff06c321d647155", size = 7033266, upload-time = "2026-09-16T09:32:05.946Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/13/f1/a2da3b55acd4ab737bf728c97edaaed5ec1d3c1236acb639dcdfa97e42c7/pypdf-6.16.2-py3-none-any.whl", hash = "sha256:c8b09a59399062fb45a1b8156c18a787a10a3dae03ac9674397a226712c94604", size = 385060, upload-time = "2026-08-23T13:50:05.349Z" }, + { url = "https://files.pythonhosted.org/packages/3c/2c/c43c03eaf630435f023f1dc61ec4a4a78951ad5530a62c71cc89bde307b7/pypdf-6.19.0-py3-none-any.whl", hash = "sha256:7e5d6e730e7dae87d560a2cee218b852f6498c8be61966f3cd02ead971e48d14", size = 395480, upload-time = "2026-09-16T09:32:04.087Z" }, ] [[package]] @@ -3115,14 +3115,14 @@ wheels = [ [[package]] name = "python-discovery" -version = "1.5.1" +version = "1.6.1" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "filelock" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/04/b7/1581a8103855c43567776aa34135e5ec3c597346c23bfd10c7eb5e0b10a4/python_discovery-1.5.1.tar.gz", hash = "sha256:e2ea8b884cd1701f386eda8cf327b87743f1dc21b7f784470799537d95635384", size = 77200, upload-time = "2026-07-31T22:06:02.48Z" } +sdist = { url = "https://files.pythonhosted.org/packages/0c/57/250bd238b966cece44328235eb85290045d059265fdaf7527a3a958123db/python_discovery-1.6.1.tar.gz", hash = "sha256:cf87d3627dfb4412437fdd5b13eae402607722998d21567993aedbc59b23c15e", size = 84338, upload-time = "2026-09-18T01:31:53.971Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/6a/07/a89b539750a159d5101c4eb9fc84e2961f65cefbd5e0b7440b284471c0b0/python_discovery-1.5.1-py3-none-any.whl", hash = "sha256:ac07f44cade589d954e9d6a1e1468539fdddd2cf676beb51da73e0f156b7c932", size = 35752, upload-time = "2026-07-31T22:06:01.116Z" }, + { url = "https://files.pythonhosted.org/packages/16/7d/e9ffbadfbf89c93848412d04594135c4ae8c1d37d9e053b9c3ed718fabc4/python_discovery-1.6.1-py3-none-any.whl", hash = "sha256:d43fcdef879fe795352bd13ccf8d185ba5a9f86f36cfcd00529f596e737442b3", size = 38664, upload-time = "2026-09-18T01:31:52.448Z" }, ] [[package]] @@ -3916,11 +3916,11 @@ wheels = [ [[package]] name = "urllib3" -version = "2.7.0" +version = "2.8.0" source = { registry = "https://pypi.org/simple" } -sdist = { url = "https://files.pythonhosted.org/packages/53/0c/06f8b233b8fd13b9e5ee11424ef85419ba0d8ba0b3138bf360be2ff56953/urllib3-2.7.0.tar.gz", hash = "sha256:231e0ec3b63ceb14667c67be60f2f2c40a518cb38b03af60abc813da26505f4c", size = 433602, upload-time = "2026-05-07T16:13:18.596Z" } +sdist = { url = "https://files.pythonhosted.org/packages/e3/05/b17359e1cefb4f909b5e40b1b90a496d987258916dbbf88e842c729f510e/urllib3-2.8.0.tar.gz", hash = "sha256:63bf2ead4c879426ebf22ef2a781eeb4aa3b4ae798a0435506f8687fd5bb9b63", size = 458972, upload-time = "2026-09-15T19:29:36.253Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/7f/3e/5db95bcf282c52709639744ca2a8b149baccf648e39c8cc87553df9eae0c/urllib3-2.7.0-py3-none-any.whl", hash = "sha256:9fb4c81ebbb1ce9531cce37674bbc6f1360472bc18ca9a553ede278ef7276897", size = 131087, upload-time = "2026-05-07T16:13:17.151Z" }, + { url = "https://files.pythonhosted.org/packages/92/9d/c4e665119135114480843e7ab388fa94d8480650450e6f8e26b70d323a4c/urllib3-2.8.0-py3-none-any.whl", hash = "sha256:0cf3cae568d36aa9576b28dfb35f11328f1cb974ca7647d9475ebb86c75ac6e3", size = 135717, upload-time = "2026-09-15T19:29:34.577Z" }, ] [[package]] @@ -3938,17 +3938,18 @@ wheels = [ [[package]] name = "virtualenv" -version = "21.7.1" +version = "21.14.6" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "distlib" }, { name = "filelock" }, + { name = "packaging" }, { name = "platformdirs" }, { name = "python-discovery" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/ea/fa/18004e5cb15541ad2a68ff219c755233b012b12d4ec8663d06a258082bec/virtualenv-21.7.1.tar.gz", hash = "sha256:d0dbfaa5483487baea28d7210ef8d24c9d1bd0f10f449eeb215568825a9b334e", size = 5525237, upload-time = "2026-07-30T15:40:36.36Z" } +sdist = { url = "https://files.pythonhosted.org/packages/92/f3/589727d02bc832750cfa00ced4315d127535a22a6d7cfcb369b2fe219a2e/virtualenv-21.14.6.tar.gz", hash = "sha256:c6b74616e64bffa9532cfa9dfd54ffbab0dc9b999df2bb8ec0d75e559b983661", size = 5467961, upload-time = "2026-10-08T14:49:57.376Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/4b/a7/ded126c19495158a05c7202b3389139839d4cf78d622d453867778e0f7a8/virtualenv-21.7.1-py3-none-any.whl", hash = "sha256:6394973f990536e34c05157179146c020284c42fe01da1dfeb0ba16c345280d9", size = 5504576, upload-time = "2026-07-30T15:40:34.512Z" }, + { url = "https://files.pythonhosted.org/packages/c9/1f/e5182163bcee149dda4eec823d88b0f51f2c0c531b12da896f666f0b80c1/virtualenv-21.14.6-py3-none-any.whl", hash = "sha256:eb1a20c5878569b69667f6a1af03458a8512b3afa7be771a2dafe92fdca05bb9", size = 5487469, upload-time = "2026-10-08T14:49:54.988Z" }, ] [[package]]