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 favor.

The diagram and score decomposition match R’s reliabilitydiag on calibrated, overconfident, squashed, heavily tied, and rare-event datasets within a tight floating-point tolerance.

Reliability Diagram¶

class calibre.evaluation.ReliabilityDiagram(prediction_values, event_probabilities, prediction_weights)[source]¶

A fitted CORP reliability diagram.

Parameters:

prediction_values (np.ndarray)
event_probabilities (np.ndarray)
prediction_weights (np.ndarray)

prediction_values¶: The distinct forecast values, ascending.

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

prediction_weights¶: Evaluation mass 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(y, x).plot(density="none")
>>> ax.get_xlabel()
'forecast probability'

calibre.corp_reliability(y_true, y_pred, *, 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:

y_true (ndarray) – Binary outcomes in {0, 1} from evaluation data not used to fit the forecaster.
y_pred (ndarray) – Forecast probabilities in [0, 1] on the same evaluation data.
sample_weight (ndarray | None) – Non-negative frequency or evaluation-mass weights. Zero-weight observations are ignored. 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(y, x)

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

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

Notes

The unweighted estimator follows Dimitriadis, Gneiting & Jordan (2021) and is pinned against their reliabilitydiag R package. sample_weight is calibre’s frequency-weight extension: integer weights give exactly the same diagram as literal row replication.

Score Decomposition¶

calibre.score_decomposition(y_true, y_pred, *, score='brier', sample_weight=None)[source]¶

Decompose a mean score into miscalibration, discrimination and uncertainty.

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

Miscalibration is what recalibration would save you, discrimination is what your forecasts buy over always predicting the base rate, and uncertainty is the difficulty of the problem, which no forecaster can change.

Parameters:

y_true (np.ndarray) – Binary outcomes in {0, 1} from evaluation data not used to fit the forecaster or calibrator.
y_pred (np.ndarray) – Forecast probabilities in [0, 1] on the same evaluation data.
score (str | Callable[[np.ndarray, np.ndarray], np.ndarray]) – Proper scoring rule: "brier" (default), "log", or a vectorized callable accepting (y_true, y_pred) and returning one finite loss per observation. The caller is responsible for ensuring a custom score is proper for binary probability forecasts.
sample_weight (np.ndarray | None) – Non-negative frequency or evaluation-mass weights. Zero-weight observations are ignored. Defaults to 1.

Returns:

mean_score, miscalibration, discrimination, and: uncertainty. Miscalibration and discrimination are non-negative for proper scoring rules, guaranteed by the optimality of the PAV solution.

Return type:

dict

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(y, x)
>>> bool(d["miscalibration"] < 0.01), bool(d["discrimination"] > 0.05)
(True, True)

The identity holds exactly:

>>> reconstructed = (
...     d["miscalibration"]
...     - d["discrimination"]
...     + d["uncertainty"]
... )
>>> bool(abs(d["mean_score"] - reconstructed) < 1e-12)
True

Warning

Evaluate held-out or out-of-fold predictions. For an isotonic-family calibrator, in-sample miscalibration is zero by construction because the calibrator and this diagnostic use the same idempotent PAV projection.

Notes

This follows Dimitriadis, Gneiting & Jordan (2021) and is pinned against their reliabilitydiag R package for the Brier score. Evaluation weights are calibre’s extension; integer weights are equivalent to row replication.

Uncertainty Bands¶

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

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

Pointwise null bands for a perfectly calibrated reliability diagram.

Forecast rows are sampled with replacement from their empirical distribution. For each sample, outcomes are drawn as y* ~ Bernoulli(y_pred) and the CORP diagram is refit. Pointwise percentiles of those refits show how far the estimated diagram can wander under the calibration hypothesis.

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

Parameters:

y_pred (ndarray) – Forecast probabilities in [0, 1]. Observed outcomes are not needed because the calibration null supplies their distribution.
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:

prediction_values (the original forecast support), lower and: upper pointwise bounds.

Return type:

dict

Notes

These are pointwise bands, not a simultaneous envelope or a global calibration test.

The method follows Dimitriadis, Gneiting & Jordan (2021), “Stable reliability diagrams for probabilistic classifiers”, and their reliabilitydiag R implementation. Only its resampling method is implemented here; the paper also derives asymptotic bands.

Each resampled fit contributes only on its observed forecast support. This also applies to singleton fits, avoiding extrapolation to forecast values absent from that resample.

calibre.confidence_bands(y_true, y_pred, *, level=0.9, n_resamples=1000, random_state=0)[source]¶

Pointwise bootstrap intervals around a CORP reliability estimate.

Forecast rows are sampled with replacement. Outcomes are then drawn from the PAV-recalibrated probabilities and a CORP diagram is refit. The returned bounds are pointwise percentiles of those refits, with the boundary correction used by the authors’ reference implementation.

Parameters:

y_true (ndarray) – Binary outcomes in {0, 1}.
y_pred (ndarray) – Forecast probabilities 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:

prediction_values (the original forecast support), lower and: upper pointwise bounds.

Return type:

dict

Warning

Evaluate held-out or out-of-fold predictions. These are pointwise bootstrap intervals, not a simultaneous envelope. Their finite-sample coverage can differ from level and should be validated for the intended data-generating setting.

Notes

This implements the resampling confidence-region method of Dimitriadis, Gneiting & Jordan (2021), “Stable reliability diagrams for probabilistic classifiers”, following their reliabilitydiag R implementation. Each resampled fit contributes only on its observed forecast support, including when that fit contains a single forecast value.

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(labels, x)
    print(
        f"{name:14s} Brier {d['mean_score']:.4f} = "
        f"MCB {d['miscalibration']:.4f} - "
        f"DSC {d['discrimination']:.4f} + UNC {d['uncertainty']:.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    "
    f"{score_decomposition(labels, in_sample)['miscalibration']:.4f}"
)
print(
    f"MCB out-of-fold  "
    f"{score_decomposition(labels, out_of_fold)['miscalibration']:.4f}"
)

The in-sample number is zero no matter how badly the model generalizes. 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).