CORP Evaluation¶

A binned reliability diagram makes the analyst pick the bins, and the picture changes with the choice. The CORP approach of Dimitriadis, Gneiting & Jordan (PNAS 2021) removes the choice: conditional event probabilities are estimated by isotonic regression via PAV, so the algorithm determines the number and position of the flat segments and there is nothing left to tune in your favour.

These numbers are pinned against R’s reliabilitydiag on five datasets (calibrated, overconfident, squashed, heavily tied, rare-event) to 1e-16 or better.

Reliability Diagram¶

class calibre.evaluation.ReliabilityDiagram(x, cep, weight)[source]¶

A fitted CORP reliability diagram.

Parameters:

x (np.ndarray)
cep (np.ndarray)
weight (np.ndarray)

x¶: The distinct forecast values, ascending.

cep¶: PAV-recalibrated conditional event probability at each forecast value.

weight¶: Number of observations carrying each forecast value.

Notes

Points where the diagram is flat are the CORP bins: the PAV algorithm chose them, so no bin count needs to be supplied and none can be tuned to flatter the forecaster.

as_function()[source]¶

Return the recalibration map as a callable.

Returns:

Piecewise-linear interpolation of: the diagram, matching the paper’s display convention. A single distinct forecast value gives a step function, since there is nothing to interpolate between.

Return type:

PiecewiseLinear or StepFunction

plot(**kwargs)[source]¶

Draw this diagram.

Convenience wrapper around calibre.plots.plot_reliability_diagram(), which documents every keyword. Needs matplotlib: pip install 'calibre[plots]'.

Parameters:: **kwargs (Any) – Passed straight through to plot_reliability_diagram() – ax, bands, density, style, diagonal, color and label.
Returns:: The axes drawn on.
Return type:: matplotlib.axes.Axes

Examples

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> import numpy as np
>>> from calibre import corp_reliability
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 1, 200)
>>> y = rng.binomial(1, x).astype(float)
>>> ax = corp_reliability(x, y).plot(density="none")
>>> ax.get_xlabel()
'forecast probability'

calibre.corp_reliability(x, y, sample_weight=None)[source]¶

Estimate conditional event probabilities by isotonic regression.

This is the CORP reliability diagram: the PAV-recalibrated forecast probabilities plotted against the original forecast values. Unlike a binned diagram it needs no bin count, because PAV determines the number and position of the flat segments itself.

Parameters:

x (ndarray) – Forecast probabilities.
y (ndarray) – Binary outcomes in {0, 1}.
sample_weight (ndarray | None) – Non-negative per-observation weights. Defaults to 1.

Returns:

The fitted diagram.

Return type:

ReliabilityDiagram

Examples

>>> import numpy as np
>>> from calibre.evaluation import corp_reliability
>>> x = np.array([0.2, 0.4, 0.6, 0.8])
>>> y = np.array([0.0, 1.0, 0.0, 1.0])
>>> diagram = corp_reliability(x, y)

The middle pair violates monotonicity, so PAV pools it to its mean:

>>> diagram.cep
array([0. , 0.5, 0.5, 1. ])

Score Decomposition¶

calibre.score_decomposition(x, y, score='brier', sample_weight=None)[source]¶

Decompose a mean score into miscalibration, discrimination and uncertainty.

Returns the CORP decomposition mean_score = MCB - DSC + UNC, where the calibrated forecasts are the PAV-recalibrated probabilities and the reference forecast is the marginal event frequency.

Read it as: MCB is what recalibration would save you, DSC is what your forecasts buy over always predicting the base rate, and UNC is the difficulty of the problem, which no forecaster can change.

Parameters:

x (ndarray) – Forecast probabilities.
y (ndarray) – Binary outcomes in {0, 1}.
score (str) – Proper scoring rule: "brier" (default) or "log".
sample_weight (ndarray | None) – Non-negative per-observation weights. Defaults to 1.

Returns:

mean_score, MCB, DSC, UNC. MCB and DSC: are non-negative, guaranteed by the optimality of the PAV solution.

Return type:

dict

Raises:

ValueError – If score is not a supported proper scoring rule.

Examples

>>> import numpy as np
>>> from calibre.evaluation import score_decomposition
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 1, 2000)
>>> y = rng.binomial(1, x).astype(float)

These forecasts are calibrated by construction, so miscalibration is small while discrimination is substantial:

>>> d = score_decomposition(x, y)
>>> bool(d["MCB"] < 0.01), bool(d["DSC"] > 0.05)
(True, True)

The identity holds exactly:

>>> bool(abs(d["mean_score"] - (d["MCB"] - d["DSC"] + d["UNC"])) < 1e-12)
True

Uncertainty Bands¶

Both are resampling-based. The paper’s asymptotic route is not implemented.

calibre.consistency_bands(x, y, level=0.9, n_resamples=1000, random_state=0)[source]¶

Bands showing how a calibrated forecaster’s diagram would scatter.

Outcomes are redrawn as y* ~ Bernoulli(x), taking the original forecasts at face value, and the diagram is refit each time. The bands therefore sit around the diagonal and answer: if these forecasts were perfectly calibrated, how far from the diagonal would the estimate wander by chance alone? An observed diagram leaving the band is the analogue of a small p-value.

Use confidence_bands() instead for an interval around the estimate.

Parameters:

x (ndarray) – Forecast probabilities.
y (ndarray) – Binary outcomes in {0, 1}. Used for the grid and for validation; the bands themselves are generated under the calibration hypothesis and do not depend on the observed outcomes.
level (float) – Nominal coverage, default 0.9.
n_resamples (int) – Number of resamples, default 1000.
random_state (int | None) – Seed. Defaults to 0 so results are reproducible.

Returns:

x (the grid), lower and upper.

Return type:

dict

Raises:

ValueError – If level is outside (0, 1) or n_resamples is below 2.

Notes

These are pointwise bands, not a simultaneous envelope. Coverage holds at each forecast value separately. It does not hold at all of them at once, and the difference is not subtle: on perfectly calibrated data the observed diagram leaves a nominal 90% band somewhere on essentially every sample. Measured over 150 replications (tests/test_monte_carlo.py):

`n`	pointwise coverage	simultaneous coverage
300	90.1%	1.3%
1200	89.6%	0.0%
4800	89.4%	0.0%

So “my curve stayed inside the band, therefore it is calibrated” is a test with a false-positive rate near one. Read the band at a forecast value you care about, or count excursions and compare that count against the nominal miss rate – do not read it as an envelope.

Resampling only. The paper also derives asymptotic bands from isotonic regression theory (a Chernoff limit for continuous forecasts), which is not implemented here; at large sample sizes this function is the expensive option rather than the unavailable one.

calibre.confidence_bands(x, y, level=0.9, n_resamples=1000, random_state=0)[source]¶

Bands around the estimated conditional event probabilities.

Outcomes are redrawn from the PAV-recalibrated probabilities rather than the original forecasts, so the bands cluster around the CORP estimate and carry the usual frequentist reading: over repeated experiments, about level of such bands contain the true conditional event probability.

Parameters:

x (ndarray) – Forecast probabilities.
y (ndarray) – Binary outcomes in {0, 1}.
level (float) – Nominal coverage, default 0.9.
n_resamples (int) – Number of resamples, default 1000.
random_state (int | None) – Seed. Defaults to 0 so results are reproducible.

Returns:

x (the grid), lower and upper.

Return type:

dict

Raises:

ValueError – If level is outside (0, 1) or n_resamples is below 2.

Notes

Pointwise, not simultaneous, exactly as for consistency_bands(); see the table there.

Coverage of the truth is below nominal on small samples. These bands are centred on the PAV-recalibrated estimate, and isotonic regression is biased at finite sample size, so the band is centred slightly off the true conditional event probability curve. Coverage of that true curve, measured against a known data-generating process over 150 replications at a nominal 90%:

`n`	coverage of the true curve
300	78.2%
1200	86.9%
4800	90.5%

The shortfall is a property of centring on an isotonic fit, not a defect in the resampling, and it vanishes as the centring bias does. Treat a 90% band on a few hundred observations as closer to an 80% one.

Resampling only; see consistency_bands().

Usage¶

Decomposing a proper score¶

mean_score = MCB - DSC + UNC holds exactly, and both MCB and DSC are non-negative by construction:

import numpy as np

from calibre import score_decomposition

rng = np.random.default_rng(0)
scores = rng.uniform(0, 1, 3000)
labels = rng.binomial(1, scores).astype(float)
overconfident = np.clip(1.6 * (scores - 0.5) + 0.5, 0, 1)

for name, x in (("honest", scores), ("overconfident", overconfident)):
    d = score_decomposition(x, labels)
    print(
        f"{name:14s} Brier {d['mean_score']:.4f} = "
        f"MCB {d['MCB']:.4f} - DSC {d['DSC']:.4f} + UNC {d['UNC']:.4f}"
    )

MCB is what recalibration would save you, DSC is what your scores buy over always predicting the base rate, and UNC is the difficulty of the problem, which no forecaster can change. A plain Brier score tells you the model got worse; this tells you which part you can fix.

Measuring honestly¶

Scoring a calibrator on the data it was fit to does not merely flatter it. For any isotonic-family calibrator it reports perfect calibration by construction, because the calibrator and the diagnostic are the same PAV projection and PAV is idempotent:

import numpy as np

from calibre import IsotonicCalibrator, cross_val_calibrate, score_decomposition

rng = np.random.default_rng(0)
scores = rng.uniform(0, 1, 1500)
labels = rng.binomial(1, scores).astype(float)

in_sample = IsotonicCalibrator().fit(scores, labels).transform(scores)
out_of_fold = cross_val_calibrate(IsotonicCalibrator(), scores, labels, cv=5)

print(f"MCB in-sample    {score_decomposition(in_sample, labels)['MCB']:.4f}")
print(f"MCB out-of-fold  {score_decomposition(out_of_fold, labels)['MCB']:.4f}")

The in-sample number is zero no matter how badly the model generalises. Use cross_val_calibrate() for any number you intend to believe.

References¶

Dimitriadis, T., Gneiting, T. & Jordan, A. I. (2021), “Stable reliability diagrams for probabilistic classifiers”, PNAS 118(8).