From 60010c991a408a8e885c46c96aa702200ce2b992 Mon Sep 17 00:00:00 2001 From: Stefan Jansen Date: Thu, 8 Oct 2026 14:53:57 -0400 Subject: [PATCH 1/2] feat(metrics): add calibration diagnostics for probability forecasts Add ml4t.diagnostic.metrics.calibration: brier_score, log_loss with clipping, brier_decomposition (Murphy reliability/resolution/uncertainty plus the Stephenson et al. within-bin terms so the identity is exact for binned continuous forecasts), reliability_table with Wilson intervals or a seeded cluster bootstrap for outcomes correlated within an event, expected/max calibration error, and calibration_slope (linear Mincer-Zarnowitz or logit Cox slope) with CR1/HC1 sandwich errors for measuring the favorite-longshot bias. Default bins are finer near 0 and 1 for prediction-market prices. Tests plant a favorite-longshot distortion and check recovery, check the decomposition identity, compare the sandwich with statsmodels, and check that the cluster bootstrap and cluster-robust errors widen for shared outcomes. Adds a Methods docs page with an executed example. Co-Authored-By: Claude Opus 5.5 --- docs/api/index.md | 9 + docs/methods/calibration.md | 188 ++++ docs/methods/index.md | 1 + docs/reference/references.md | 16 + mkdocs.yml | 1 + src/ml4t/diagnostic/metrics/AGENTS.md | 3 + src/ml4t/diagnostic/metrics/__init__.py | 22 + src/ml4t/diagnostic/metrics/calibration.py | 807 ++++++++++++++++++ tests/metrics/test_calibration.py | 421 +++++++++ .../metrics/test_public_metrics_namespace.py | 7 + 10 files changed, 1475 insertions(+) create mode 100644 docs/methods/calibration.md create mode 100644 src/ml4t/diagnostic/metrics/calibration.py create mode 100644 tests/metrics/test_calibration.py 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__)) From 6572fde0b0851367a30d030bb2e4c43537f44d9f Mon Sep 17 00:00:00 2001 From: Stefan Jansen Date: Thu, 8 Oct 2026 15:07:07 -0400 Subject: [PATCH 2/2] fix(deps): lock patched pypdf, urllib3 and virtualenv The PyPI advisory audit in the security workflow reported 19 known vulnerabilities in pypdf 6.16.2, urllib3 2.7.0 and virtualenv 21.7.1, failing every Audit job. Lock pypdf 6.19.0, urllib3 2.8.0 and virtualenv 21.14.6 (python-discovery follows virtualenv); pip-audit now reports no known vulnerabilities. Co-Authored-By: Claude Opus 5.5 --- uv.lock | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) 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]]