Skip to content

Metrics and tests

Data these snippets assume

Blocks on this page use a fixed set of names for held-out data instead of re-deriving it each time. To run them, define the names first:

# docs: no-run — the definitions the snippets on this page assume
import numpy as np
from probcal import BetaCalibrator, Masterscale, make_pd_portfolio
from probcal.monitor import CalibrationMonitor

cal_set, new_set = make_pd_portfolio(n=3000, random_state=0), make_pd_portfolio(n=1000, random_state=1)
s_cal, y_cal = cal_set.scores, cal_set.y          # held-out calibration scores and outcomes
w_cal = np.ones_like(y_cal)                       # uniform sample weights
s_new = new_set.scores                            # scores of new obligors
ms = Masterscale.from_edges([0.01, 0.05], names=["G1", "G2", "G3"])   # the masterscale
p_cal = BetaCalibrator().fit(s_cal, y_cal).predict_proba(s_cal)         # calibrated PDs
grades = ms.assign(p_cal)                                              # rating labels
segments = np.array(["seg-a", "seg-b", "seg-c"])[np.arange(len(s_cal)) % 3]  # segment labels
model = ...   # any object with predict_proba(X); the docs use a stub that reads X[:, 0]
mon = CalibrationMonitor(alpha=0.05)              # with a few batches of a calibrated forecast applied

Every block also names, in its first comment line, which of these it uses.

Measuring calibration is harder than fixing it. The quantity of interest, \( \Pr(Y = 1 \mid \hat{p}) \), is a conditional expectation that no finite sample reveals directly, so every metric estimates it through some smoothing device, and every smoothing device imports bias, sensitivity, or both. This chapter walks the full catalog implemented in probcal.metrics, states each estimator's formula and known pathologies, and ends with the table that answers the operational question: which of these may I select a calibrator on, and which are for reporting only? The short answer, argued in detail below: select on log loss (default) or Brier; report the ECE family and ICI; never select on ECE or Hosmer–Lemeshow.

All metrics share the signature metric(y, p, *, sample_weight=None, **kw), and evaluate(y, p) assembles everything into a MetricReport with seeded bootstrap percentile confidence intervals.

Sample weights

Every sample_weight in probcal.metrics and probcal.curves must be strictly positive and finite, and is read as a relative frequency weight:

  • Point estimates (scores, binned and smooth calibration errors, fitted intercept and slope) are weighted means. They do not change when all weights are multiplied by a constant, and unit weights reproduce the unweighted numbers exactly.
  • Wherever a sample size enters (a test statistic's variance, a finite-sample bias correction, a null envelope, a likelihood-ratio statistic), the variance is the one the weighted statistic actually has when each y is Bernoulli(p): Var(sum(w * (y - p))) = sum(w**2 * p * (1 - p)). spiegelhalter_z and the ecce_curve null envelope use that variance directly; hosmer_lemeshow uses it per risk group, in the classical form sum_g(w**2) * pt_g * (1 - pt_g) with pt_g = sum_g(w**2 p) / sum_g(w**2); ece_debiased and murphy_decomposition(bias_corrected=True) use each bin's effective count sum_b(w)**2 / sum_b(w**2); and the likelihood-ratio tests of calibration_test and calibration_belt are taken on the weighted log-likelihood and referred to a second-order Rao–Scott corrected chi-square, with sandwich covariances for the belt's bands (so CalibrationTestResult.df is fractional when the weights vary). Every one of these is invariant to multiplying all weights by a constant, and unit weights give the unweighted statistic exactly.

One count-form function reads weights literally instead: pluto_tasche_from_arrays sums weights as obligor and default counts (an integer weight equals row duplication). hl_e_test raises each likelihood-ratio factor to min(1, w / mean(w)) within its grade, matching CalibrationMonitor, so it too is invariant to rescaling the weights. The per-grade binomial and Jeffreys tests use raw integer counts and ignore weights.

Because the test statistics use the exact weighted variances above, Hosmer–Lemeshow and the calibration LR test do not change when every weight is multiplied by a constant c. A single global effective sample size such as Kish's sum(w)**2 / sum(w**2) would be scale-invariant too, but it misstates the variance whenever the weights vary with p. For weights that are independent of p the two agree closely; for weights correlated with p only the exact weighted variances give correctly sized tests.

Proper scoring rules

The log loss and Brier score (Brier, 1950) were defined in Why calibration:

\[ \mathrm{LL} = -\frac{1}{n}\sum_i \bigl[y_i \ln p_i + (1-y_i)\ln(1-p_i)\bigr], \qquad \mathrm{BS} = \frac{1}{n}\sum_i (p_i - y_i)^2 . \]

Both are strictly proper: their expectation is uniquely minimized by the true conditional probability, so no calibration map can improve them by lying. Their sample versions are unbiased estimators of the expected loss, which is the property none of the direct calibration metrics below share, and the reason they anchor selection. Log loss penalizes tail overconfidence harshly (a confident wrong prediction costs unboundedly); the Brier score is bounded and correspondingly gentler. For a portfolio-level orientation the Brier skill score

\[ \mathrm{BSS} = 1 - \frac{\mathrm{BS}}{\mathrm{BS}_{\text{ref}}}, \qquad \mathrm{BS}_{\text{ref}} = \bar{y}(1 - \bar{y}), \]

references the climatology forecast \( p \equiv \bar{y} \): positive values beat the base rate, and on a 3% portfolio the reference is small, so seemingly modest Brier differences are large skill differences.

Murphy decomposition. The binned estimator of Murphy's (1973) partition,

\[ \mathrm{BS} \approx \underbrace{\frac{1}{n}\sum_b n_b\,(\bar{p}_b - \bar{y}_b)^2}_{\text{reliability}} \;-\; \underbrace{\frac{1}{n}\sum_b n_b\,(\bar{y}_b - \bar{y})^2}_{\text{resolution}} \;+\; \underbrace{\bar{y}(1-\bar{y})}_{\text{uncertainty}}, \]

with \( \bar p_b, \bar y_b \) the mean prediction and event rate in bin \( b \), makes the calibration–sharpness trade visible in two numbers. It inherits every bias of its binning: the plug-in reliability term is biased upward and resolution downward, exactly the effect Bröcker (2009) formalized and Ferro and Fricker (2012) corrected. probcal implements the corrected variant alongside the naive one and documents the binning dependence rather than hiding it.

The Murphy diagram. murphy_curve computes a complementary, binning-free view of the same Brier score: the elementary score of the "act if \( p > \theta \)" decision rule,

\[ S_\theta(p, y) = \theta \cdot \mathbb{1}\{p > \theta,\, y = 0\} + (1 - \theta) \cdot \mathbb{1}\{p \le \theta,\, y = 1\}, \]

whose weighted mean over \( \theta \in [0, 1] \) traces the diagram; doubling its integral recovers the Brier score exactly (Ehm, Gneiting, Jordan and Krüger, 2016), since a single observation's continuous integral is \( p^2/2 \) (\( y = 0 \)) or \( (1-p)^2/2 \) (\( y = 1 \)), exactly half the Brier contribution. Plotting \( S_\theta(A) - S_\theta(B) \) for two forecasts (plots.plot_murphy(..., diff=True)) shows where along the decision spectrum one beats the other, rather than collapsing the comparison to a single Brier difference that a cancellation across thresholds can mask. murphy_curve defaults thresholds=513 (numpy.linspace(0, 1, 513), the package's dense-grid convention), evaluated in \( O(n \log n + T \log n) \) by sorting \( p \) once; the discrete identity converges to the exact one at rate roughly \( 1/n \). \( S_\theta \) is piecewise linear between consecutive unique \( p \) values but jumps exactly there, so a trapezoid over the default grid recovers Brier to about \( 10^{-3} \) on typical portfolios, tightening as \( n \) grows. Isotonic (PAV) recalibration never increases \( S_\theta \) at any threshold, so a raw forecast plotted against its own PAV fit diagnoses the value of recalibration pointwise across the whole decision range instead of in one scalar.

The analogous calibration–refinement split of the log loss replaces squared gaps with Kullback–Leibler terms: with \( c(p) = \Pr(Y=1 \mid \hat p = p) \) estimated by a recalibration curve, calibration is the mean divergence between \( \mathrm{Bernoulli}(c(p)) \) and \( \mathrm{Bernoulli}(p) \), refinement the mean entropy of \( \mathrm{Bernoulli}(c(p)) \). The split is only as good as the plug-in estimate of \( c \); the estimator choice is recorded in the changelog when implemented.

Binned estimators

The expected calibration error family discretizes the conditional expectation with \( B \) bins:

\[ \mathrm{ECE} = \sum_{b=1}^{B} \frac{n_b}{n}\,\bigl|\bar{p}_b - \bar{y}_b\bigr|, \]

with variants: strategy="mass" (equal-count bins, the recommended default) or "width"; norm="l2" squares the gaps; norm="max" takes the worst bin, which is the maximum calibration error (MCE). Two pathologies are structural, not incidental. First, binning sensitivity: ECE is a function of \( B \) and the bin edges, and rankings of models can flip under a different, equally defensible binning. Second, finite-sample bias: within each bin the absolute difference of two noisy means is biased upward, so a perfectly calibrated model has positive expected ECE, and the bias grows with \( B \) and shrinks with portfolio size slowly. ece_debiased applies the bias correction in the spirit of Bröcker (2009) and Ferro and Fricker (2012); ece_sweep implements the monotonic-sweep calibration error of Roelofs et al. (2022), a principled, data-driven resolution choice that markedly reduces bias. It follows the paper's rule (its eq. 8 and Algorithm 1): sweep \( B = 2, 3, \ldots \) equal-mass bins and stop at the first bin count whose bin event rates are not monotone; the chosen \( B^* \) is the last count before that violation, so every smaller count is monotone too. The paper sweeps up to \( B = n \); probcal stops at max_bins=100 as a cost cap. ece_sweep(rule="largest") instead takes the largest monotone \( B \) anywhere in 2..100, skipping over non-monotone counts. On synthetic continuous data the two rules choose a different \( B \) in about half of the cases (the first-violation rule usually the smaller one), and the values differ by at most ~0.01 absolute at n=500 and ~0.003 at n=10⁴.

The Hosmer–Lemeshow test (Hosmer and Lemeshow, 1980) groups observations into \( g \) risk deciles (hosmer_lemeshow(..., n_bins=10)) and forms

\[ C = \sum_{b=1}^{g} \frac{(O_b - E_b)^2}{E_b\,\bigl(1 - E_b/n_b\bigr)} \;\sim\; \chi^2_{g - \text{ddof}}, \]

with \( O_b \) observed and \( E_b \) expected events per group. The degrees of freedom depend on where \( p \) came from. For a logistic model fitted on the same data, the classical \( g-2 \) applies (ddof=2). When \( p \) is fixed before the outcomes are seen, which is the usual PD validation setting with a model or calibrator fitted on other data, the per-group residuals are asymptotically independent and \( C \sim \chi^2_g \) (Hosmer, Lemeshow and Sturdivant, 2013, on external validation). That is the default, ddof=0. Against \( g-2 \), out-of-sample \( p \) rejects about 12% of the time at a nominal 5% with \( g=10 \). Under \( \chi^2_g \) it rejects about 5%. With heavy-tailed weights and few effective events per group, the chi-square approximation is loose: with lognormal(0, 1) weights, it rejects about 9% at n = 5000 and PD ≈ 3%, and 6% at n = 50 000. It is the historical workhorse of clinical model validation and it carries the same two diseases in sharper form: the statistic depends on an essentially arbitrary grouping (changing \( g \), or the tie handling at decile boundaries, changes the p-value), and its power scales with \( n \) so that on large portfolios it rejects calibration defects of no practical consequence, while on small ones it detects almost nothing. probcal ships it because validators expect it, marks it report-only, and never lets the selector see it. When no more than ddof non-empty groups remain, the reference distribution has no degrees of freedom left, so the result reports df=0 and p_value=nan rather than a p-value on a floored df=1.

Binning-free estimators

Four estimators avoid the binning choice altogether.

Smooth ECE (Błasiok and Nakkiran, 2024) replaces hard bins with kernel smoothing: the residuals \( y_i - p_i \) are smoothed with a reflected Gaussian kernel (probcal applies it on the logit scale) and the calibration error is read from the smoothed curve, with the bandwidth chosen by the paper's self-consistency principle: the reported error is the fixed point where the measurement scale matches the error magnitude. The result is a continuous quantity, insensitive to reparametrization, with none of ECE's edge artifacts; any implementation simplification is recorded in the changelog.

ECCE, the empirical cumulative calibration error (Arrieta-Ibarra, Gujral, Tannen, Tygert and Xu, 2022), sorts observations by \( p \) and tracks the cumulative deviation \( C_k = \sum_{i \le k} (y_{(i)} - p_{(i)}) \). Under calibration this walk is a martingale hovering near zero; systematic over- or under-prediction makes it drift. The Kolmogorov-style maximum \( \max_k |C_k|/n \) and the mean absolute deviation summarize the drift, and the plot of \( C_k \) against sorted \( p \) localizes where the miscalibration lives without any smoothing parameter at all; ecce_curve and plot_ecce render it (Visualization).

Ties. Sorting by \( p \) leaves the order inside a block of tied predictions arbitrary, and reading the walk inside such a block would make the statistic depend on input row order. probcal reads the walk only at the end of each tied block (where it is order-free), takes the \( x \) coordinate as the cumulative weight fraction, and averages the mean deviation over blocks weighted by block weight. Without ties and with unit weights this is the textbook statistic, bit for bit; with ties ecce_max differs from the textbook walk only if its maximum falls inside a block, while ecce_mean differs (about +8.5% on scores rounded to two decimals), and EcceCurve has one point per distinct prediction. ecce, ecce_curve and evaluate share the same walk.

ICI and its quantiles (Austin and Steyerberg, 2019). Fit a LOESS smoother \( \hat{c}(p) \) of outcome on prediction (Austin and Steyerberg, 2014, established the graphical practice) and average the absolute distance to the diagonal:

\[ \mathrm{ICI} = \frac{1}{n} \sum_i \bigl|\hat{c}(p_i) - p_i\bigr| , \]

with E50, E90 and Emax the median, 90th percentile, and maximum of the same distances. The family is smooth, interpretable in probability units, and inherits only the mild LOESS-bandwidth dependence (frac=0.75 by default, stated openly).

Spiegelhalter's z (Spiegelhalter, 1986) is the classical unbiasedness test built directly on the Brier score. Its numerator \( \sum_i (y_i - p_i)(1 - 2p_i) \) has expectation zero under calibration, and standardizing by its variance under the null,

\[ z = \frac{\sum_i (y_i - p_i)(1 - 2p_i)} {\sqrt{\sum_i (1 - 2p_i)^2\, p_i (1 - p_i)}} , \]

gives an asymptotically standard normal statistic with a two-sided p-value. No binning, no smoothing; the trade is that it aggregates over the whole range and can miss compensating regional errors. When every prediction is exactly 0.5 the null variance is zero and the statistic is undefined: z and p_value are NaN (and the guardrail reports spiegelhalter_ok=False).

Kernel calibration error and tests (SKCE)

The squared kernel calibration error (Widmann, Lindsten and Zachariah, 2019) embeds the residual measure in a reproducing-kernel Hilbert space: the population SKCE is zero exactly when the model is calibrated, for any universal kernel. In probcal's binary specialization (the paper's identity-matrix kernel construction with predictions represented as \( (1-p, p) \)), the kernel term reduces to

\[ h_{ij} = 2\, \tilde{k}(s_i, s_j)\, (y_i - p_i)(y_j - p_j) , \]

where the factor 2 keeps values comparable with the paper's framework. The residuals \( y_i - p_i \) always stay on the probability scale; only the kernel input \( s \) may be logit-transformed (scale="logit", the low-PD option). Three consistent estimators (the paper's Table 1):

Estimator Definition Properties
"biased" \( n^{-2} \sum_{i,j} h_{ij} \) (diagonal included) V-statistic; squared RKHS norm, always ≥ 0; biased upward
"uq" (default) \( (n(n-1))^{-1} \sum_{i \ne j} h_{ij} \) unbiased; may be negative
"ul" \( \lfloor n/2 \rfloor^{-1} \sum_i h_{(2i-1),(2i)} \) over seeded disjoint pairs unbiased; O(n); higher variance

The defaults follow the paper's own experiments: a Laplacian kernel \( \exp(-|d|/\mathrm{bw}) \) with the median-heuristic bandwidth (Gretton et al., 2012), implemented deterministically (an evenly strided subsample of at most 4096 points above \( n = 4096 \), a mean-distance fallback when heavy ties drive the median to zero, and a refusal with instructions when all scores are identical). A Gaussian kernel is available.

skce_test turns the estimate into a one-sided test of H0: calibrated. The default method="bootstrap" uses the quadratic statistic with the Arcones–Giné (1992) centered resampling, the construction the paper itself states in its Appendix G rather than the wild bootstrap other implementations substitute, at O(n_boot · n²) cost. method="asymptotic" uses the linear estimator with a normal approximation at O(n), the practical choice for \( n \gtrsim 20\,000 \); the trade, stated openly, is power. A single random pairing can miss slope-type miscalibration (residual means that change sign across the score range) that the bootstrap test rejects on the same data, the paper's documented power gap. Both report p_value_bound, the distribution-free bound \( \min\!\bigl(1, \exp(-\lfloor n/2 \rfloor\, t^2 / 8)\bigr) \): valid without any asymptotics but loose, so treat it as a worst-case check and decide with p_value.

Against Spiegelhalter's z the contrast is scope: the z tests one global moment condition and can miss compensating regional errors, while the SKCE is sensitive to any deviation the kernel can resolve at its bandwidth. Neither skce nor skce_test accepts sample_weight: the U-statistic theory behind unbiasedness, the bootstrap, and the bounds is stated for unweighted i.i.d. samples, and probcal refuses to improvise weighted inference the source does not cover. Kumar, Sarawagi and Jain's (2018) MMCE is a special case of the SKCE (the paper's Example I.1), so it is not implemented separately.

The recalibration-regression framework

The most decision-relevant diagnostics come from Cox's (1958) idea of regressing the outcome on the prediction. Three quantities, all fitted by the shared IRLS core:

The calibration intercept fits \( \operatorname{logit} \Pr(Y=1) = \alpha + \operatorname{logit}(p) \) with the slope fixed at 1 (an offset-term logistic regression): \( \alpha \) is calibration-in-the-large in log-odds. On a PD portfolio, \( \alpha = -0.3 \) says the model overestimates portfolio risk by a factor \( e^{0.3} \approx 1.35 \) in odds.

The calibration slope fits \( \operatorname{logit} \Pr(Y=1) = \alpha + \beta\, \operatorname{logit}(p) \) and reads \( \beta \): values below 1 mean predictions are too spread out (the signature of overfitting), and values above 1 mean underfitting. These are the same parameters a Platt calibrator would fit as repairs, here estimated as diagnoses.

The calibration test is the likelihood-ratio test of \( (\alpha, \beta) = (0, 1) \) jointly, on 2 degrees of freedom: the Cox-framed "weak calibration" test, in the lineage running through Miller, Hui and Tierney (1991). Its χ² p-value comes from probcal._math.chi2_sf, an upper-tail evaluation that keeps the runtime numpy-only and stays accurate deep in the tail (every χ² and normal p-value in the package, including Hosmer–Lemeshow, the belt, Spiegelhalter and SKCE, is computed this way, so tiny p-values are reported as such instead of rounding to exactly 0).

Guardrails. calibration_guardrails(y, p) condenses the framework into three flags used across the package and printed in every selection report: slope within \( [0.9, 1.1] \), intercept within \( \pm 0.1 \), Spiegelhalter p-value above 0.05. The thresholds are conventions, not theorems; they encode "no deviation a validator would flag" and are documented as such.

Per-grade backtesting

Credit-risk validation operates on rating grades, not on continuous scores. Given grade assignments and per-grade PDs, binomial_grade_test computes for each grade the exact binomial tail probability of observing at least the realized number of defaults under the grade's PD (the incomplete-beta representation via probcal._math.betainc keeps it exact at any \( n \)), alongside the normal approximation, in the traffic-light style summary supervisors expect (BCBS, 2005). jeffreys_grade_test implements the ECB's preferred formulation (ECB, 2019): the posterior for the grade's true default rate under the Jeffreys prior is \( \mathrm{Beta}(k + \tfrac12,\; n - k + \tfrac12) \), and the reported p-value is the posterior probability that the true rate lies at or below the assigned PD. The reading is one-sided and conservative by design, so a small value flags a grade whose PD is likely understated. The documentation says so explicitly, because two-sided misreadings of the Jeffreys test are a recurring validation error. Both results carry 90% display intervals (ci_low/ci_high) for plot_grade_backtest (Visualization); the intervals are for reading, the traffic lights carry the verdict.

Uncertainty: the bootstrap protocol

A metric without an uncertainty statement invites overreading, and calibration metrics on percent-level event rates are noisy in ways intuition underestimates. evaluate(y, p) therefore attaches confidence intervals to every scalar it reports, by the case-resampling bootstrap: draw \( n \) observations with replacement from the evaluation pairs, recompute the metric, repeat n_boot=1000 times with a seeded generator, and report the 2.5th and 97.5th percentiles of the bootstrap distribution. The percentile method is chosen because it is simple to implement and easy to state: it makes no normality assumption, which matters for bounded and skewed statistics like ECE near zero, and it is reproducible bit for bit given the seed.

Stratification is the default. Each replicate resamples the negative and positive classes separately (case resampling within strata, the pROC-style default), so every replicate reproduces the observed class counts exactly. This conditions the CI on the observed class balance: it excludes the additional variance a plain i.i.d. bootstrap picks up from the event count itself fluctuating replicate to replicate, which on low-event-rate data can be the dominant source of resampling noise. Excluding that source of variance can narrow the interval relative to i.i.d. resampling on exactly the rare-event, small-\( n \) data where it matters most. That is the opposite of the intuition that stratifying always tightens or always widens a CI: the direction depends on which variance source dominates. evaluate(..., stratify=False) restores plain i.i.d. resampling, redrawing a degenerate (single-class) replicate up to 100 times before raising rather than silently substituting anything. Reusing the point estimate as a zero-variance replicate whenever an i.i.d. draw comes back single-class would narrow the i.i.d. path artificially rather than report the sampling variance honestly. Neither the stratified default nor the redraw rule makes CIs uniformly wider or narrower; both make the reported interval mean what it claims to measure.

Bootstrap intervals for biased estimators still center on the biased value. A bootstrap CI around plain ECE quantifies its variance, not its bias, so the interval can exclude zero for a perfectly calibrated model. The report pairs ECE with its debiased variant precisely so this artifact is visible rather than misread. Bootstrap-heavy computations carry the slow pytest marker and a fixed default seed, per the package's reproducibility conventions.

Grouping (by=). evaluate(..., by=labels) runs the exact same pooled call above on the full data plus one independent call per group, returning a GroupedMetricReport instead of a plain MetricReport. Groups are formed from the raw label values (so 1 and "1" are different groups) and displayed as str(value), ordered by display name; two distinct values that print the same raise ValueError. Group i (in that order) uses random_state + 1000 * i rather than reusing random_state for every group. The offset is fixed and label-independent, so reproducibility does not depend on how many groups exist or what they are named, and no two groups' bootstrap draws can coincide by construction. This is side-by-side reporting, not a test: no comparison across groups is computed, and no multiple-comparison correction is applied, because none is implied by returning several independent reports. Formal group-conditional calibration testing is future work; see docs/guide/groups.md.

Weighted quantiles. e50, e90, and the reliability_summary stats box compute their quantile step with probcal._math.weighted_quantile whenever sample_weight is given: the continuous weighted form of numpy's default (linear) rule. Each sorted observation sits at the midpoint of its own weight, those midpoints are mapped affinely onto [0, 1] as the plotting positions, and the quantile interpolates linearly between them. With equal weights the positions are exactly numpy's, so the result is np.quantile bit for bit; the positions are continuous in the weights, so nearly equal weights give nearly the unweighted answer, and rescaling the weights changes nothing.

Reading a report

A MetricReport is designed to be read in a fixed order. Start with log loss and Brier against their pre-calibration values: did the repair help at all, and is the improvement larger than the bootstrap intervals overlap? Then the guardrail triplet of slope, intercept and Spiegelhalter, which localizes any remaining defect to spread, level, or neither. Then the descriptive family (debiased ECE, smooth ECE, ICI, ECCE), read as a cross-check: these should broadly agree, and when they do not, the disagreement itself is diagnostic (a large MCE with small ICI means one bad region, not global miscalibration; a large ECCE maximum with small mean means a localized drift). Per-grade tests come last, because they answer a different question: not "is the map good" but "which grades would a supervisor flag". The visualization chapter pairs each layer of this reading with a plot.

What to select on: the table

Metric Proper Binning-sensitive Finite-sample bias Formal test Selection use
Log loss strictly no unbiased no default criterion
Brier score strictly no unbiased no alternative criterion
Brier skill score derived no mild (ratio) no report
Murphy / LL decompositions n/a yes corrected variant available no report
ECE (mass/width, MCE) no yes upward no never
Debiased ECE no yes reduced no report
ECE sweep (Roelofs) no reduced reduced no optional, with care
Smooth ECE no no low no optional, with care
ECCE no no low max-statistic report
ICI / E50 / E90 / Emax no no (LOESS frac) low no optional
Spiegelhalter z n/a no n/a yes never (it is a test)
SKCE (skce, skce_test) n/a no (bandwidth) uq/ul unbiased yes never (it is a test)
Hosmer–Lemeshow no yes n/a yes never
Calibration intercept/slope n/a no n/a yes (LR/Wald) guardrails
Binomial / Jeffreys per grade n/a grades fixed n/a yes never (backtest)

The logic behind the verdicts compresses to one principle. Selection is optimization, and optimizing a biased, binning-dependent, non-proper quantity invites the optimizer to exploit the estimator rather than improve the calibration. A calibrator can win an ECE contest by emitting values that straddle bin edges favorably, and win a Hosmer–Lemeshow contest by blurring predictions until the test loses power. Strictly proper scores close that loophole by construction. The selector therefore defaults to out-of-fold log loss, accepts Brier, ICI, smooth ECE and ECE-sweep as deliberate alternatives, refuses plain ECE and Hosmer–Lemeshow entirely, and prints the guardrail flags next to whatever criterion was used.

In probcal

# s_cal, y_cal, ms: held-out calibration scores, outcomes, the masterscale
from probcal.metrics import (
    brier_score, calibration_guardrails, calibration_slope, ece, ece_debiased,
    evaluate, ici, jeffreys_grade_test, log_loss, skce, skce_test, smooth_ece,
    spiegelhalter_z,
)

print(log_loss(y_cal, s_cal), brier_score(y_cal, s_cal))  # proper: safe to select on
print(ece(y_cal, s_cal), ece_debiased(y_cal, s_cal))      # report-only; note the bias
print(smooth_ece(y_cal, s_cal), ici(y_cal, s_cal))        # binning-free
print(calibration_slope(y_cal, s_cal), spiegelhalter_z(y_cal, s_cal))
print(skce(y_cal, s_cal), skce_test(y_cal, s_cal).p_value)  # kernel calibration error + test
print(calibration_guardrails(y_cal, s_cal))               # the three-flag summary

report = evaluate(y_cal, s_cal, n_boot=100, random_state=42)   # everything + bootstrap CIs
print(report)

print(jeffreys_grade_test(y_cal, s_cal, ms))              # ECB-style backtest, best to worst

Computational cost

Most of the catalog is O(n) or O(n log n) per call: the proper scores, ECCE, Spiegelhalter's z, and the recalibration-regression framework are single linear passes; the binned ECE family sorts or bins in O(n log n). Two estimators smooth rather than bin, which costs more, and both take an anchoring parameter that brings their cost down without changing what they measure.

The ICI family (ici, e50, e90, emax, and the reliability_summary stats box) fits a LOESS smoother \( \hat{c}(p) \). Refitting it at every one of the \( n \) observations, each fit itself scanning an O(n)-window, is effectively O(n²) at portfolio scale. grid_size (default 512) fits the smoother at that many equal-mass anchors spanning the prediction range and linearly interpolates the rest, the same device R's stats::lowess uses via its delta parameter. Windows and bandwidths are computed against the full data, so this changes how many points get an exact fit, not what the fit means; measured drift on make_pd_portfolio(n=5000) is |Δici| ≈ 1.3e-6, far below bootstrap CI width. grid_size=None gives the exact per-point fit at that O(n²) cost. On this host, ici at n=50,000 takes 1.2s with the default grid against 192.2s exact, and loess(grid_size=512) fits n=1,000,000 points in under 30s.

smooth_ece solves a bandwidth fixed point by bisection; evaluated exactly, each step builds a kernel matrix against every residual. bins (default 8192) pre-aggregates the weighted residual measure onto equal-width bins over the logit range once, up front (O(n)); each bisection step then evaluates that binned measure directly on its own lattice by truncated-Gaussian convolution, independent of n. A small-bandwidth guard retries once on an adaptively refined binning (bins <- ceil(range / (sigma/8))) whenever the found bandwidth would be under-resolved by the current bins, then falls back to the exact per-observation computation only if that refinement is infeasible (above 2^20 bins) or still under-resolved, so accuracy never degrades silently. The binned path evaluates on its own lattice rather than the exact path's 257-point grid, which would alias against the bin lattice. The lattice path engages for every call with a non-degenerate logit range, whatever n; bins=None or a degenerate range runs the exact computation. For n <= bins the lattice value may differ from the exact grid at the ~1e-4 level on typical portfolios (measured ≤ 2.4e-4 on make_pd_portfolio); on wide clipped-logit-range data the gap can be much larger, because there the exact path's fixed 257-point grid under-resolves small-sigma kernels and the lattice value (≥ 8 samples per sigma) is the better one.

evaluate's cost is dominated by the bootstrap, not any single metric: every point estimate in the requested catalog is recomputed n_boot times (default 1000; n_boot=0 returns point estimates only, with NaN interval bounds). Per replicate, scores, regression metrics and ECCE are O(n); the binned ECEs are O(n log n), with ece/ece_debiased/mce sharing one 15-bin pass; ece_sweep is O(n log n) per candidate bin count and stops at the first non-monotone count (typically fewer than 15 candidates); the ICI family shares one LOESS fit at O(grid_size · frac · n); smooth_ece bins once in O(n) and then costs O(bins · taps) per bisection step, where taps is the truncated-Gaussian kernel width (at most ~161 taps), independent of n. metrics= restricts the catalog to the names actually needed.

Each replicate is sorted by prediction once (np.argsort(..., kind="stable")), and the LOESS fit and the ECCE walk reuse that order. The LOESS anchors are fitted in cache-sized blocks of whole windows rather than one Python iteration each (_math._loess_fit_sorted_vec, the only LOESS engine; it matches a scalar per-point loop to ≤ 4.3×10⁻¹⁵ on fitted values and is tested against it). Reported point estimates come off the unsorted path; the bootstrap replicates' reordered sums move percentile CI bounds only in their last bits.

The vectorized anchor fit differs from a scalar loop by ≤ 2.3×10⁻¹⁶ relative on well-conditioned windows only. If a window is rank-deficient (every non-zero tricube weight sitting on one distinct p, which needs the far half of the window to lie at exactly the bandwidth), the local-linear determinant is pure cancellation (~10⁻²³ rather than 0), and an ulp-level weight difference can put two evaluation orders on opposite sides of the abs(det) < _FPMIN guard, giving values that differ by O(1). On such a window the swy / sw branch (the weighted mean, what a rank-deficient local linear fit degenerates to) is the well-defined answer; the other branch divides by cancellation noise. Since anchors are data quantiles, this has not been observed to reach a reported value, and the corner is pinned by tests/test_math.py.

Measured on the dev host at n=10⁴, one full-catalog replicate costs about 28 ms; ece_sweep with its early stop costs 1.8 ms per replicate (about 8.9 ms without it), and the binned ECE family costs 0.34 ms. evaluate(n=10⁴, n_boot=1000) takes about 28 s. The ICI family's LOESS fit remains the largest single cost, so metrics= excluding it (ici/e50/e90/emax) is the biggest lever. For n above roughly 10⁶, reduce n_boot, pass a metrics= subset, or both.

docs/scripts/benchmarks.py measures single calls on demand. Measured on the dev host:

n call wall time (s)
10,000 ece(d.y, d.scores) 0.004
10,000 ece_sweep(d.y, d.scores) 0.002
10,000 ecce(d.y, d.scores) 0.001
10,000 ici(d.y, d.scores) 0.021
10,000 smooth_ece(d.y, d.scores) 0.001
10,000 evaluate(d.y, d.scores, n_boot=100) 2.8
100,000 ece(d.y, d.scores) 0.006
100,000 ece_sweep(d.y, d.scores) 0.058
100,000 ecce(d.y, d.scores) 0.010
100,000 ici(d.y, d.scores) 0.193
100,000 smooth_ece(d.y, d.scores) 0.003
100,000 evaluate(d.y, d.scores, n_boot=100) 25.8

The public ece_sweep and ici run the same code as the bootstrap path (early-stopping sweep, vectorized LOESS), so the single-call rows above are the per-replicate costs.

References

  • Arcones, M. A., Giné, E. (1992). "On the bootstrap of U and V statistics." Annals of Statistics 20(2), 655–674.
  • Arrieta-Ibarra, I., Gujral, P., Tannen, J., Tygert, M., Xu, C. (2022). "Metrics of Calibration for Probabilistic Predictions." Journal of Machine Learning Research 23(351), 1–54.
  • Austin, P. C., Steyerberg, E. W. (2014). "Graphical assessment of internal and external calibration of logistic regression models by using loess smoothers." Statistics in Medicine 33(3), 517–535.
  • Austin, P. C., Steyerberg, E. W. (2019). "The Integrated Calibration Index (ICI) and related metrics for quantifying the calibration of logistic regression models." Statistics in Medicine 38(21), 4051–4065.
  • BCBS (2005). Studies on the Validation of Internal Rating Systems. Working Paper No. 14, revised version, May 2005. Bank for International Settlements.
  • Błasiok, J., Nakkiran, P. (2024). "Smooth ECE: Principled Reliability Diagrams via Kernel Smoothing." ICLR.
  • Brier, G. W. (1950). "Verification of forecasts expressed in terms of probability." Monthly Weather Review 78(1), 1–3.
  • Bröcker, J. (2009). "Reliability, sufficiency, and the decomposition of proper scores." Quarterly Journal of the Royal Meteorological Society 135(643), 1512–1519.
  • Cox, D. R. (1958). "Two further applications of a model for binary regression." Biometrika 45, 562–565.
  • ECB (2019). Instructions for reporting the validation results of internal models: IRB Pillar I models for credit risk. European Central Bank Banking Supervision, February 2019.
  • Ehm, W., Gneiting, T., Jordan, A., Krüger, F. (2016). "Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings." Journal of the Royal Statistical Society: Series B 78(3), 505–562.
  • Ferro, C. A. T., Fricker, T. E. (2012). "A bias-corrected decomposition of the Brier score." Quarterly Journal of the Royal Meteorological Society 138(668), 1954–1960.
  • Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., Smola, A. (2012). "A Kernel Two-Sample Test." Journal of Machine Learning Research 13, 723–773.
  • Hosmer, D. W., Lemeshow, S. (1980). "Goodness of fit tests for the multiple logistic regression model." Communications in Statistics - Theory and Methods 9(10), 1043–1069.
  • Hosmer, D. W., Lemeshow, S., Sturdivant, R. X. (2013). Applied Logistic Regression, 3rd ed. Wiley.
  • Kumar, A., Sarawagi, S., Jain, U. (2018). "Trainable Calibration Measures for Neural Networks from Kernel Mean Embeddings." ICML, PMLR 80, 2805–2814.
  • Miller, M. E., Hui, S. L., Tierney, W. M. (1991). "Validation techniques for logistic regression models." Statistics in Medicine 10(8), 1213–1226.
  • Murphy, A. H. (1973). "A New Vector Partition of the Probability Score." Journal of Applied Meteorology 12(4), 595–600.
  • Roelofs, R., Cain, N., Shlens, J., Mozer, M. C. (2022). "Mitigating Bias in Calibration Error Estimation." AISTATS, PMLR 151, 4036–4054.
  • Spiegelhalter, D. J. (1986). "Probabilistic prediction in patient management and clinical trials." Statistics in Medicine 5(5), 421–433.
  • Widmann, D., Lindsten, F., Zachariah, D. (2019). "Calibration tests in multi-class classification: A unifying framework." NeurIPS 32.