A worked calibrator comparison

No calibrator is best for every population. This tutorial compares several methods on three known distortions and keeps calibration fitting separate from evaluation. It is a worked example, not a portable speed benchmark or a universal ranking.

The comparison uses held-out Brier score as the selection criterion because it is a proper score. CORP miscalibration (MCB), discrimination (DSC), smECE, and the number of distinct predictions explain why the proper score changed.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split

from calibre import (
    CenteredIsotonicCalibrator,
    IsotonicCalibrator,
    NearlyIsotonicCalibrator,
    RelaxedPAVACalibrator,
    SplineCalibrator,
    brier_score,
    score_decomposition,
    smooth_calibration_error,
)

RANDOM_STATE = 42
random_generator = np.random.default_rng(RANDOM_STATE)

Generate known distortions

The population event probability is sampled first. Outcomes come from that probability, while forecasts apply a known transformation on the log-odds scale:

• overconfident forecasts multiply the true log odds by 1.8;

• underconfident forecasts multiply them by 0.55;

• shifted forecasts add 0.8 to the true log odds.

These are monotone distortions. A separate nonmonotone design would be needed to evaluate methods intended to permit ranking violations.

def probability(log_odds):
    """Map finite log odds to probabilities."""
    return 1.0 / (1.0 + np.exp(-log_odds))


def make_scenario(scale=1.0, shift=0.0, n_observations=3000):
    """Draw outcomes and forecasts under a known log-odds distortion."""
    true_probability = random_generator.beta(2.0, 2.0, n_observations)
    outcomes = random_generator.binomial(1, true_probability).astype(float)
    true_log_odds = np.log(true_probability / (1.0 - true_probability))
    forecasts = probability(scale * true_log_odds + shift)
    return forecasts, outcomes


scenarios = {
    "overconfident": make_scenario(scale=1.8),
    "underconfident": make_scenario(scale=0.55),
    "shifted": make_scenario(shift=0.8),
}

Define comparable candidates

The candidates below are general-purpose calibrators. Relaxed PAVA needs an explicit scientific bound, so the example states one. CDI isotonic is not included because it requires application-specific decision thresholds on the forecast scale; dropping arbitrary thresholds into a generic benchmark would not evaluate its intended use.

def candidate_calibrators():
    """Return fresh calibrators for one scenario."""
    return {
        "isotonic": IsotonicCalibrator(),
        "centered isotonic": CenteredIsotonicCalibrator(),
        "nearly isotonic": NearlyIsotonicCalibrator(),
        "spline": SplineCalibrator(),
        "relaxed PAVA (-0.01)": RelaxedPAVACalibrator(min_increment=-0.01),
    }

Fit on calibration data and evaluate elsewhere

Every candidate sees the same calibration rows and the same held-out rows. Errors are allowed to stop the documentation build; silently assigning an infinite score would hide a broken calibrator.

def evaluation_row(scenario, method, outcomes, predictions):
    """Summarize one held-out prediction vector."""
    decomposition = score_decomposition(outcomes, predictions)
    return {
        "scenario": scenario,
        "method": method,
        "brier": brier_score(outcomes, predictions),
        "miscalibration": decomposition["miscalibration"],
        "discrimination": decomposition["discrimination"],
        "smece": smooth_calibration_error(outcomes, predictions),
        "distinct_predictions": np.unique(predictions).size,
    }


rows = []
for scenario_name, (forecasts, outcomes) in scenarios.items():
    (
        calibration_forecasts,
        test_forecasts,
        calibration_outcomes,
        test_outcomes,
    ) = train_test_split(
        forecasts,
        outcomes,
        test_size=0.5,
        stratify=outcomes,
        random_state=RANDOM_STATE,
    )

    rows.append(
        evaluation_row(
            scenario_name,
            "uncalibrated",
            test_outcomes,
            test_forecasts,
        )
    )

    for method_name, calibrator in candidate_calibrators().items():
        calibrator.fit(calibration_forecasts, calibration_outcomes)
        calibrated_predictions = calibrator.transform(test_forecasts)
        rows.append(
            evaluation_row(
                scenario_name,
                method_name,
                test_outcomes,
                calibrated_predictions,
            )
        )

results = pd.DataFrame(rows)
results.assign(
    brier=results["brier"].round(4),
    miscalibration=results["miscalibration"].round(4),
    discrimination=results["discrimination"].round(4),
    smece=results["smece"].round(4),
)

	scenario	method	brier	miscalibration	discrimination	smece	distinct_predictions
0	overconfident	uncalibrated	0.2024	0.0100	0.0576	0.0702	1500
1	overconfident	isotonic	0.1975	0.0020	0.0545	0.0253	56
2	overconfident	centered isotonic	0.1973	0.0049	0.0576	0.0285	1475
3	overconfident	nearly isotonic	0.1976	0.0020	0.0545	0.0253	57
4	overconfident	spline	0.1972	0.0048	0.0576	0.0278	1496
5	overconfident	relaxed PAVA (-0.01)	0.2054	0.0047	0.0492	0.0261	1458
6	underconfident	uncalibrated	0.2068	0.0112	0.0543	0.0652	1500
7	underconfident	isotonic	0.2012	0.0024	0.0511	0.0258	36
8	underconfident	centered isotonic	0.2005	0.0049	0.0543	0.0282	1439
9	underconfident	nearly isotonic	0.2012	0.0024	0.0511	0.0258	37
10	underconfident	spline	0.2003	0.0047	0.0543	0.0295	1500
11	underconfident	relaxed PAVA (-0.01)	0.2092	0.0056	0.0463	0.0380	1428
12	shifted	uncalibrated	0.2163	0.0290	0.0624	0.1504	1500
13	shifted	isotonic	0.1918	0.0025	0.0606	0.0281	57
14	shifted	centered isotonic	0.1913	0.0040	0.0624	0.0284	1447
15	shifted	nearly isotonic	0.1918	0.0026	0.0606	0.0281	60
16	shifted	spline	0.1920	0.0047	0.0624	0.0371	1499
17	shifted	relaxed PAVA (-0.01)	0.2043	0.0065	0.0520	0.0452	1420

Read the comparison

Lower Brier is the decision rule. MCB and smECE diagnose calibration, DSC measures useful separation, and distinct predictions describe resolution. Neither a low calibration error nor a large number of distinct predictions is sufficient by itself.

for scenario_name in scenarios:
    scenario_results = results[results["scenario"] == scenario_name]
    ordered = scenario_results.sort_values("brier")
    print(f"\n{scenario_name}")
    for row in ordered.itertuples():
        print(
            f"  {row.method:22s} Brier {row.brier:.4f}  "
            f"MCB {row.miscalibration:.4f}  DSC {row.discrimination:.4f}  "
            f"distinct {row.distinct_predictions}"
        )

overconfident
  spline                 Brier 0.1972  MCB 0.0048  DSC 0.0576  distinct 1496
  centered isotonic      Brier 0.1973  MCB 0.0049  DSC 0.0576  distinct 1475
  isotonic               Brier 0.1975  MCB 0.0020  DSC 0.0545  distinct 56
  nearly isotonic        Brier 0.1976  MCB 0.0020  DSC 0.0545  distinct 57
  uncalibrated           Brier 0.2024  MCB 0.0100  DSC 0.0576  distinct 1500
  relaxed PAVA (-0.01)   Brier 0.2054  MCB 0.0047  DSC 0.0492  distinct 1458

underconfident
  spline                 Brier 0.2003  MCB 0.0047  DSC 0.0543  distinct 1500
  centered isotonic      Brier 0.2005  MCB 0.0049  DSC 0.0543  distinct 1439
  nearly isotonic        Brier 0.2012  MCB 0.0024  DSC 0.0511  distinct 37
  isotonic               Brier 0.2012  MCB 0.0024  DSC 0.0511  distinct 36
  uncalibrated           Brier 0.2068  MCB 0.0112  DSC 0.0543  distinct 1500
  relaxed PAVA (-0.01)   Brier 0.2092  MCB 0.0056  DSC 0.0463  distinct 1428

shifted
  centered isotonic      Brier 0.1913  MCB 0.0040  DSC 0.0624  distinct 1447
  isotonic               Brier 0.1918  MCB 0.0025  DSC 0.0606  distinct 57
  nearly isotonic        Brier 0.1918  MCB 0.0026  DSC 0.0606  distinct 60
  spline                 Brier 0.1920  MCB 0.0047  DSC 0.0624  distinct 1499
  relaxed PAVA (-0.01)   Brier 0.2043  MCB 0.0065  DSC 0.0520  distinct 1420
  uncalibrated           Brier 0.2163  MCB 0.0290  DSC 0.0624  distinct 1500

The lowest value in each block describes this simulated split only. Repeating the simulation with new populations, or using repeated validation splits on real data, is necessary before making a stable method choice.

figure, axes = plt.subplots(1, 3, figsize=(16, 4.5))

for axis, column, title, log_scale in (
    (axes[0], "brier", "Held-out Brier score", False),
    (axes[1], "miscalibration", "CORP miscalibration", False),
    (axes[2], "distinct_predictions", "Distinct held-out predictions", True),
):
    pivoted = results.pivot(index="scenario", columns="method", values=column)
    pivoted.plot.bar(ax=axis, legend=False)
    axis.set_title(title)
    axis.set_xlabel("")
    axis.tick_params(axis="x", rotation=0)
    if log_scale:
        axis.set_yscale("log")

handles, labels = axes[0].get_legend_handles_labels()
figure.legend(handles, labels, loc="outside lower center", ncol=3)
figure.tight_layout(rect=(0, 0.15, 1, 1))
plt.show()

Choosing a method in an application

Use the same sequence on real data:

1. Fit each candidate on calibration observations the classifier did not train on.

2. Compare candidates on separate validation observations with log loss or Brier score.

3. Inspect held-out calibration error, discrimination, and resolution to understand the choice.

4. After choosing, evaluate once on a final test set that played no role in fitting or selection.

Isotonic and centered isotonic fit once with no tuning search. The automatic nearly isotonic and spline fits cross-validate their regularization choices, so they generally cost more. Runtime depends on data size, candidate grids, software versions, and hardware; measure it on the intended workload rather than copying a timing from this tutorial.