mixle.inference.scoring module

Proper scoring rules for probabilistic forecasts.

A scoring rule S(F, y) measures how well a predictive distribution F matched the value y that actually occurred; it is proper when the forecaster minimises the expected score by reporting their true belief (and strictly proper when that minimiser is unique). These give a single, fair currency for comparing probabilistic forecasts and interval methods at matched coverage – a 90% interval that is wide everywhere and a 90% interval that is tight where it can be both cover 90% of the time, but the second wins on the interval score.

Lower is better for every score here (they are penalties / losses):

  • log_score() – the log loss -log p(y) (local, strictly proper).

  • brier_score() / brier_decomposition() – squared-error score for categorical forecasts, with the reliability/resolution/uncertainty decomposition.

  • crps_ensemble() / crps_gaussian() – the Continuous Ranked Probability Score, the CDF-space analogue of absolute error; reduces to absolute error for a point forecast.

  • interval_score() (a.k.a. winkler_score()) – the proper score for a central interval forecast: width plus an out-of-bounds penalty scaled by the miscoverage level.

  • pinball_loss() – the check loss whose expectation is minimised by the true quantile.

  • energy_score() – the multivariate generalisation of CRPS.

  • skill_score() – a score expressed as fractional improvement over a reference.

All functions are pure NumPy/SciPy and operate on plain arrays; by default they return the mean score over observations, with mean=False returning the per-observation vector for paired comparisons (see mixle.inference.resampling() for CIs on score differences).

log_score(prob, *, mean=True)[source]

Logarithmic score (log loss) -log p(y) of the probability assigned to the outcome.

The log score is local (it depends only on the density at the realised value) and strictly proper. Pass the predictive probability/density evaluated at the observed outcome.

Parameters:

  • prob (ndarray) – array of predictive probabilities (or densities) at the realised outcomes. Must be finite and non-negative – a probability or density can never be negative, NaN, or infinite. Zero is a legitimate (if poor) forecast, and is clipped away from zero before the log so it yields a large but finite penalty rather than inf.

  • mean (bool) – if True (default) return the mean log loss; otherwise the per-observation vector.

Returns:

Mean log loss (float) or the per-observation array.

Raises:

ValueError – if any entry of prob is negative, NaN, or infinite. Clipping a negative value the same way as the legitimate zero-probability edge case would let an invalid input (e.g. -1) silently masquerade as a valid, if poor, forecast; and an infinite input would silently become the score -inf – which, this module being lower-is-better, looks like the single best score achievable, letting a broken upstream density silently win any comparison instead of being rejected.

Return type:

ndarray | float

brier_score(prob, outcome, *, mean=True)[source]

Brier score (mean squared error of probabilistic classification).

For binary forecasts pass 1-D prob (predicted probability of the positive class) and 0/1 outcome. For K-class forecasts pass prob shaped (n, K) and outcome either as integer class labels (length n) or as a one-hot (n, K) matrix; the multiclass score sums the squared error across classes, so it ranges in [0, 2].

Parameters:

  • prob (ndarray) – (n,) positive-class probabilities, or (n, K) class probabilities. Must be finite and lie in [0, 1] – these are genuine probabilities, not densities.

  • outcome (ndarray) – (n,) 0/1 or integer labels, or (n, K) one-hot. Must be finite, and each entry a legitimate label for prob’s shape: 0/1 against a 1-D binary prob, an integer class index in [0, K) against an (n, K) prob, or 0/1 against an (n, K) one-hot matrix.

  • mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean Brier score (float) or the per-observation array.

Raises:

ValueError – if prob is not finite or has entries outside [0, 1], or if outcome is not finite or has entries outside its valid domain. Without these checks an invalid prob/outcome pair can cancel out into a deceptively “perfect” score of 0.0 (e.g. brier_score([-1, 2], [-1, 2])), and an out-of-range integer label silently corrupts an unrelated one-hot slot instead of raising.

Return type:

ndarray | float

brier_decomposition(prob, outcome, *, bins=10)[source]

Murphy’s reliability–resolution–uncertainty decomposition of the (binary) Brier score.

Bins the forecast probabilities into bins equal-width bins on [0, 1] and computes

Brier = reliability - resolution + uncertainty,

where reliability (lower is better) is the mean squared gap between forecast probability and observed frequency within a bin, resolution (higher is better) is how far bin frequencies move from the base rate, and uncertainty is the base-rate variance p_bar (1 - p_bar) (a property of the data, not the forecaster).

Parameters:

  • prob (ndarray) – (n,) predicted probabilities of the positive class. Must be finite and lie in [0, 1].

  • outcome (ndarray) – (n,) 0/1 outcomes. Must be finite and each entry exactly 0 or 1.

  • bins (int) – number of equal-width probability bins.

Returns:

{'reliability', 'resolution', 'uncertainty', 'brier', 'brier_binned_reconstruction'}. brier is the ACTUAL Brier score mean((prob - outcome)^2) (audit SCORING-3: the key previously held the binned reconstruction, which differs from the real score by the within-bin variance term). brier_binned_reconstruction = reliability - resolution + uncertainty is the decomposition’s own sum. The components are PLUG-IN estimates whose finite-sample bias grows like bins/n (audit SCORING-2): with sparse bins, reliability is inflated by within-bin binomial noise and resolution by base-rate noise – read them against n/bins, and do not chase small differences at small samples.

Raises:

ValueError – if prob is not finite or has entries outside [0, 1], or if outcome is not finite or has entries outside {0, 1}. Without these checks an invalid input pair can silently produce a mathematically impossible negative “brier” value instead of raising (e.g. brier_decomposition([-1, 2], [-1, 2]) returns brier == -2.0).

Return type:

dict[str, float]

crps_ensemble(forecasts, y, *, fair=True, mean=True)[source]

Continuous Ranked Probability Score from a finite predictive ensemble (sample).

Uses the energy form CRPS = E|X - y| - 1/2 E|X - X'| estimated from the ensemble draws. The CRPS generalises absolute error to distributional forecasts (a point forecast recovers |x-y|) and is reported in the same units as y.

Parameters:

  • forecasts (ndarray) – predictive draws. Either (n, m) (m draws for each of n observations) or (m,) for a single observation. Ragged ensembles are not supported – pad or call per observation.

  • y (ndarray) – (n,) realised values (or a scalar for the (m,) case).

  • fair (bool) – if True (the DEFAULT – audit SCORING-1) use the unbiased 1/(m(m-1)) spread estimator (the “fair” CRPS), which is proper for finite ensembles: measured at m = 5, the expected fair score is minimized by reporting the true predictive (dispersion factor 1.000), while the 1/m^2 estimator rewards UNDER-dispersing to a factor of 0.700 – an improper score that trains and selects forecasters toward overconfidence. fair=False keeps the plain 1/m^2 estimator for comparability with sources that use it; treat those numbers as biased-low in ensemble spread.

  • mean (bool) – if True return the mean CRPS; otherwise the per-observation vector.

Returns:

Mean CRPS (float) or the per-observation array.

Return type:

ndarray | float

crps_gaussian(mu, sigma, y, *, mean=True)[source]

Closed-form CRPS for a Gaussian predictive distribution N(mu, sigma^2).

CRPS = sigma * [ z (2 Phi(z) - 1) + 2 phi(z) - 1/sqrt(pi) ] with z = (y - mu) / sigma (Gneiting & Raftery 2007). Exact, so it is the right reference when the forecast is Gaussian.

Parameters:

  • mu (ndarray) – predictive means, broadcastable with y.

  • sigma (ndarray) – predictive standard deviations (> 0), broadcastable with y.

  • y (ndarray) – realised values.

  • mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean CRPS (float) or the per-observation array.

Return type:

ndarray | float

interval_score(lower, upper, y, alpha, *, mean=True)[source]

Winkler interval score for a central (1 - alpha) prediction interval.

IS = (u - l) + (2/alpha)(l - y) 1[y < l] + (2/alpha)(y - u) 1[y > u] (Gneiting & Raftery 2007). It rewards narrow intervals but penalises misses by 2/alpha times the shortfall, so at matched nominal coverage the tighter, better-placed interval scores lower. This is the proper way to rank interval methods (conformal, quantile-regression, GP credible bands) head to head.

Parameters:

  • lower (ndarray) – (n,) interval lower endpoints.

  • upper (ndarray) – (n,) interval upper endpoints.

  • y (ndarray) – (n,) realised values.

  • alpha (float) – miscoverage level (e.g. 0.1 for a 90% interval), in (0, 1).

  • mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean interval score (float) or the per-observation array.

Return type:

ndarray | float

winkler_score(lower, upper, y, alpha, *, mean=True)[source]

Alias for interval_score() (the score is due to Winkler 1972).

Parameters:
Return type:

ndarray | float

pinball_loss(pred, y, tau, *, mean=True)[source]

Pinball (quantile / check) loss for a quantile forecast.

rho_tau(u) = max(tau * u, (tau - 1) * u) with u = y - pred; its expectation is minimised by the true tau-quantile, so it is the proper score for quantile forecasts and the loss that quantile regression and conformalised quantiles optimise.

Parameters:

  • pred (ndarray) – predicted tau-quantiles, shape (n,) (single level) or (n, q) (several levels, one per column of tau).

  • y (ndarray) – (n,) realised values.

  • tau (float | ndarray) – quantile level(s) in (0, 1) – a scalar for (n,) pred, or a length-q array matching the columns of a (n, q) pred.

  • mean (bool) – if True return the mean over all entries; otherwise the per-observation (or per-observation-per-level) array.

Returns:

Mean pinball loss (float) or the per-observation array.

Return type:

ndarray | float

energy_score(forecasts, y, *, fair=True, mean=True)[source]

Energy score: the multivariate generalisation of CRPS for vector-valued forecasts.

ES = E||X - y|| - 1/2 E||X - X'|| with Euclidean norms; for scalar y it equals crps_ensemble().

Parameters:

  • forecasts (ndarray) – predictive draws shaped (n, m, d) (m d-vectors per observation) or (m, d) for a single observation.

  • y (ndarray) – (n, d) realised vectors (or (d,) for the single-observation case).

  • fair (bool) – if True (the DEFAULT – audit SCORING-1, same finite-ensemble impropriety as crps_ensemble()) use the unbiased 1/(m(m-1)) spread estimator; False keeps the biased 1/m^2 form for comparability, rewarding under-dispersion.

  • mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean energy score (float) or the per-observation array.

Return type:

ndarray | float

skill_score(score, reference, *, perfect=0.0)[source]

Skill score: fractional improvement of score over a reference score.

skill = (reference - score) / (reference - perfect) – 1 means the forecast hit the perfect score, 0 means it only matched the reference (e.g. climatology), and negative means it did worse than the reference. Works for any negatively-oriented score (lower is better) such as those above.

Parameters:

  • score (float) – the forecast’s score (lower is better).

  • reference (float) – the baseline/reference forecast’s score.

  • perfect (float) – the score of a perfect forecast (0.0 for the rules in this module).

Returns:

The skill score (float); nan if the reference already equals the perfect score.

Return type:

float