Probability Calibration

A probability model is calibrated when cases assigned probability 0.8 are positive 80% of the time in the population, and the same agreement holds at its other probability levels. A finite sample only estimates this agreement. Calibration differs from ranking: a model can rank cases well while systematically overstating their probabilities. The classification-metrics article uses probabilities to compare decision costs; here we examine whether those probabilities are reliable.

Reading a reliability curve

A reliability curve groups predictions into bins and plots each bin’s mean predicted probability on the horizontal axis and observed positive fraction on the vertical axis. Points above the diagonal indicate under-prediction, and points below indicate over-prediction. For example, a bin averaging 0.25 with 12 positives among 100 cases has an observed frequency of 0.12 and a gap of −0.13. Bin frequencies fluctuate with sample size; a point off the diagonal is an estimate, not proof of a population gap. Report bin counts when assessing how precise those estimates are.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.naive_bayes import GaussianNB
from sklearn.calibration import calibration_curve
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(0)
X0, y = make_classification(n_samples=30000, n_features=8, n_informative=5,
                            n_redundant=0, class_sep=0.8, random_state=0)
# 20 noisy copies violate naive Bayes conditional independence
X = np.hstack([X0, X0[:, [0]] + rng.normal(0, 0.35, (30000, 20))])
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0)

p = GaussianNB().fit(Xtr, ytr).predict_proba(Xte)[:, 1]
observed, predicted = calibration_curve(yte, p, n_bins=10, strategy="uniform")
counts = np.bincount(np.searchsorted(np.linspace(0, 1, 11)[1:-1], p), minlength=10)
for count, o, q in zip(counts[counts > 0], observed, predicted):
    print(f"n {count:5d}   predicted {q:.4f}   observed {o:.4f}   gap {o - q:+.4f}")
# n  2840   predicted 0.0301   observed 0.0359   gap +0.0059
# n  1100   predicted 0.1472   observed 0.1082   gap -0.0390
# n   894   predicted 0.2499   observed 0.1230   gap -0.1269
# n   827   predicted 0.3493   observed 0.1995   gap -0.1498
# n   896   predicted 0.4515   observed 0.2980   gap -0.1535
# n  1054   predicted 0.5521   observed 0.4858   gap -0.0663
# n  1367   predicted 0.6540   observed 0.6452   gap -0.0088
# n  1962   predicted 0.7518   observed 0.8048   gap +0.0530
# n  2163   predicted 0.8517   observed 0.9011   gap +0.0493
# n  1897   predicted 0.9386   observed 0.9805   gap +0.0419

Gaussian naive Bayes models features as conditionally independent given the class. The 20 noisy copies are dependent on the original feature, so multiplying their likelihood contributions treats related measurements as separate evidence. They are not 21 identical likelihood factors, and this experiment does not isolate duplication as the sole cause of the curve. In the displayed bins, predicted probabilities around 0.15–0.55 exceed observed frequencies by as much as 0.154, while bins above 0.75 underpredict by about 0.04–0.05. A signed average gap could cancel opposing errors; an absolute-gap summary would not cancel them across bins.

Brier, log loss, and what they weigh

Both are strictly proper scoring rules: their population expected loss is uniquely minimized by the true conditional probability. That property specifies a target; finite data, a restricted model, or regularization can still produce miscalibrated fitted probabilities. Brier and log loss assess probability quality, but neither isolates calibration from the model’s ability to distinguish cases. Use reliability diagnostics alongside them, and evaluate decision costs separately when those costs define the application.

Log loss penalizes extreme errors more heavily. For an actual positive, a report of 0.01 gives Brier loss 0.9801 and log loss about 4.605; at 0.000001 these become about 1 and 13.816. Their numerical scales differ, so ratios between the two metrics are not meaningful. In exercise 2, log loss increases more than ninefold between two calibration procedures, while Brier increases by about 15%. Neither percentage alone specifies an application’s cost.

Platt scaling and isotonic regression

A calibrator is another fitted model that maps a base score to a probability. Fit it on predictions for cases excluded from base-model fitting, then evaluate it on separate data. Sigmoid calibration fits \(r(s)=1/(1+\exp[-(as+b)])\); isotonic regression fits a non-decreasing mapping, allowing more flexibility. Sigmoid calibration is strictly increasing when \(a>0\); its fitted slope is not constrained to be positive in this implementation. A decreasing map would reverse the ranking. Isotonic fits pooled levels at observed scores; scikit-learn interpolates between fitted thresholds for new scores.

from sklearn.calibration import CalibratedClassifierCV
from sklearn.ensemble import RandomForestClassifier
from sklearn.frozen import FrozenEstimator          # sklearn >= 1.6; replaces cv="prefit"
from sklearn.metrics import brier_score_loss, log_loss, roc_auc_score

Xcal, Xte2, ycal, yte2 = train_test_split(Xte, yte, test_size=0.5, random_state=0)
rf = RandomForestClassifier(n_estimators=300, random_state=0, n_jobs=-1).fit(Xtr, ytr)

p_raw = rf.predict_proba(Xte2)[:, 1]
print(f"{'raw':9s} AUC {roc_auc_score(yte2, p_raw):.6f}"
      f"  Brier {brier_score_loss(yte2, p_raw):.4f}  logloss {log_loss(yte2, p_raw):.4f}")
for method in ("sigmoid", "isotonic"):
    cal = CalibratedClassifierCV(FrozenEstimator(rf), method=method).fit(Xcal, ycal)
    p = cal.predict_proba(Xte2)[:, 1]
    print(f"{method:9s} AUC {roc_auc_score(yte2, p):.6f}"
          f"  Brier {brier_score_loss(yte2, p):.4f}  logloss {log_loss(yte2, p):.4f}")
# raw       AUC 0.974098  Brier 0.0690  logloss 0.2500
# sigmoid   AUC 0.974098  Brier 0.0590  logloss 0.1996
# isotonic  AUC 0.973380  Brier 0.0587  logloss 0.2031

On the separate evaluation subset, sigmoid reduces log loss from 0.2500 to 0.1996, about 20%, and preserves the reported AUC. Its increasing fitted map preserves ranking mathematically, although numerical saturation can create ties. Isotonic has the slightly lower Brier score here, 0.0587 against 0.0590, but higher log loss, 0.2031 against 0.1996. Its AUC falls from 0.974098 to 0.973380 as scores are pooled into ties. Such ties can raise or lower AUC depending on which positive–negative pairs are merged; the decrease is an observation from this run.

The code requires scikit-learn 1.6 or later for FrozenEstimator, which keeps the already-fitted forest fixed while the calibrator is fitted. Here the original split supplies 15,000 base-training rows; the remaining rows are divided into 7,500 calibration and 7,500 evaluation rows. The examples reuse data for demonstration, so they are not a sequence of independent final tests. In an actual development workflow, choose the calibration method without repeatedly consulting the final test set. Passing an unfitted estimator with integer cv supports cross-validated calibration; with ensemble=True, it fits a base-model/calibrator pair on complementary portions in each fold and averages their predictions. The split must still respect groups and time.

ECE is not a single number

Expected calibration error (ECE) here is the sample-weighted mean absolute bin gap: \(\sum_b(n_b/n)|\bar y_b-\bar p_b|\). A bin containing 20% of the cases with a gap of 0.1 contributes 0.02. Uniform bins divide the probability range into equal widths; quantile bins aim for similar case counts. Repeated scores can make quantile edges coincide, leaving fewer nonempty bins than requested.

def ece(y, p, n_bins, strategy="uniform"):
    edges = (np.linspace(0, 1, n_bins + 1) if strategy == "uniform"
             else np.quantile(p, np.linspace(0, 1, n_bins + 1)))
    idx = np.clip(np.searchsorted(edges, p, side="right") - 1, 0, n_bins - 1)
    return sum((idx == b).mean() * abs(y[idx == b].mean() - p[idx == b].mean())
               for b in range(n_bins) if (idx == b).sum())

p_nb = GaussianNB().fit(Xtr, ytr).predict_proba(Xte)[:, 1]
for nb in (5, 10, 20, 50, 100, 500):
    print(f"bins {nb:4d}   ECE uniform {ece(yte, p_nb, nb):.4f}"
          f"   ECE quantile {ece(yte, p_nb, nb, 'q'):.4f}")
# bins    5   ECE uniform 0.0499   ECE quantile 0.0530
# bins   10   ECE uniform 0.0538   ECE quantile 0.0534
# bins   20   ECE uniform 0.0564   ECE quantile 0.0561
# bins   50   ECE uniform 0.0573   ECE quantile 0.0566
# bins  100   ECE uniform 0.0586   ECE quantile 0.0580
# bins  500   ECE uniform 0.0692   ECE quantile 0.0712

On these predictions, uniform-bin ECE rises from 0.0499 with five bins to 0.0692 with 500, so the five-bin figure is about 28% lower. This change does not identify a single cause. Coarse bins can hide variation within a bin, while fine bins have fewer observations and noisier positive fractions. Finite-sample variation and binning can bias ECE estimates; even calibrated probabilities can yield nonzero empirical ECE. The opposing gaps below 0.6 and above 0.75 are mostly in different coarse bins, so that pattern alone does not explain the increase. With bins that are not nested, changing their number can also change their boundaries; ECE need not rise monotonically in general.

Report the bin count, edge strategy, and evaluation sample size with ECE. Inspect the reliability curve and the number of observations supporting each region. Brier and log loss avoid binning, but they measure overall probability quality and therefore answer a different question. A useful evaluation can include both kinds of evidence.

Two models, one dataset, different calibration

The next example fits gradient boosting and a random forest on the same credit-default training data and uses identical probability-band boundaries. The models differ in both ranking and probability estimation, and their bands contain different accounts. The fairlearn dataset fetch needs internet access on its first call.

import numpy as np
from fairlearn.datasets import fetch_credit_card
from sklearn.model_selection import train_test_split
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.metrics import roc_auc_score, brier_score_loss

cc = fetch_credit_card(as_frame=True)
X = cc.data.astype(float).to_numpy()
y = cc.target.astype(int).to_numpy()
Xtr, Xte, ytr, yte = train_test_split(
    X, y, test_size=0.3, random_state=0, stratify=y)

hgb = HistGradientBoostingClassifier(random_state=0).fit(Xtr, ytr).predict_proba(Xte)[:, 1]
rf = RandomForestClassifier(n_estimators=200, random_state=0,
                            n_jobs=-1).fit(Xtr, ytr).predict_proba(Xte)[:, 1]

for name, q in (("HistGradientBoosting", hgb), ("RandomForest", rf)):
    print(f"{name:22s} AUC {roc_auc_score(yte, q):.4f}  Brier {brier_score_loss(yte, q):.4f}")
    for lo, hi in [(0.0, 0.1), (0.1, 0.3), (0.3, 0.5), (0.5, 0.7), (0.7, 1.01)]:
        k = (q >= lo) & (q < hi)
        if k.sum() > 30:
            print(f"   band {lo:.1f}-{hi:.1f}  n {k.sum():5,}  predicted {q[k].mean():.4f}"
                  f"  observed {yte[k].mean():.4f}  gap {yte[k].mean() - q[k].mean():+.4f}")
# HistGradientBoosting   AUC 0.7765  Brier 0.1347
#    band 0.0-0.1  n 2,578  predicted 0.0656  observed 0.0671  gap +0.0015
#    band 0.1-0.3  n 4,431  predicted 0.1681  observed 0.1690  gap +0.0009
#    band 0.3-0.5  n   938  predicted 0.3775  observed 0.3785  gap +0.0009
#    band 0.5-0.7  n   587  predicted 0.6113  observed 0.6201  gap +0.0088
#    band 0.7-1.0  n   466  predicted 0.7588  observed 0.7511  gap -0.0078
# RandomForest           AUC 0.7627  Brier 0.1392
#    band 0.0-0.1  n 2,427  predicted 0.0572  observed 0.0746  gap +0.0174
#    band 0.1-0.3  n 4,314  predicted 0.1763  observed 0.1627  gap -0.0136
#    band 0.3-0.5  n 1,117  predicted 0.3861  observed 0.3491  gap -0.0370
#    band 0.5-0.7  n   707  predicted 0.5944  observed 0.5573  gap -0.0371
#    band 0.7-1.0  n   435  predicted 0.7907  observed 0.7448  gap -0.0459

The boosted model’s five observed band gaps are small, with maximum absolute gap 0.0088. This is evidence at the chosen bin resolution on one test sample, not calibration everywhere. Broad bands can hide within-band patterns or subgroup differences. Log-loss training targets conditional probabilities, but does not guarantee a calibrated fitted model.

The forest underpredicts in the lowest band by 0.0174 and overpredicts in the other four bands by 0.0136–0.0459. In its top band, the average report is 0.7907 and the observed default fraction is 0.7448. These are empirical default frequencies, not classification accuracies. Whether that band drives a decision depends on the policy and its costs.

The boosted model has better AUC and Brier scores in this run. Their gaps are on different scales and cannot be compared as measures of how much ranking or calibration deteriorated. Similar overall AUCs do not establish similar precision among the top few percent; evaluate that subset directly if it is the operating point. For expected monetary loss, case probabilities must be combined with the relevant loss amounts. Positive and negative probability errors may offset or accumulate depending on their signs and those amounts; these band averages do not determine the portfolio error.

Exercises

1. Which calibrator, and why. Compare Platt scaling and isotonic regression across calibration-set sizes, under two different kinds of miscalibration.

You should get: different average Brier rankings under the two distortions. Consider the calibration family, sample size, and information retained in the score.

import numpy as np
from sklearn.isotonic import IsotonicRegression
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import brier_score_loss

def logit(p):
    p = np.clip(p, 1e-6, 1 - 1e-6)
    return np.log(p / (1 - p))

def run(distort, label):
    print(f"\n{label}")
    print(f"{'n_cal':>7} {'raw':>8} {'sigmoid':>9} {'isotonic':>9}   winner")
    for n_cal in (200, 1000, 5000, 20000):
        s_, i_, r_ = [], [], []
        for rep in range(15):
            g = np.random.default_rng(rep)
            p_true = g.beta(2, 2, n_cal); y = (g.random(n_cal) < p_true).astype(int)
            score = distort(p_true)
            pt = g.beta(2, 2, 40000); yt = (g.random(40000) < pt).astype(int)
            st = distort(pt)
            r_.append(brier_score_loss(yt, np.clip(st, 1e-6, 1 - 1e-6)))
            lr = LogisticRegression().fit(logit(score).reshape(-1, 1), y)
            s_.append(brier_score_loss(yt, lr.predict_proba(logit(st).reshape(-1, 1))[:, 1]))
            iso = IsotonicRegression(out_of_bounds="clip").fit(score, y)
            i_.append(brier_score_loss(yt, iso.predict(st)))
        s, i, r = np.mean(s_), np.mean(i_), np.mean(r_)
        print(f"{n_cal:7d} {r:8.4f} {s:9.4f} {i:9.4f}   {'sigmoid' if s < i else 'isotonic'}")

run(lambda p: 1 / (1 + np.exp(-(2.5 * logit(p) + 0.8))), "A: sigmoidal distortion")
run(lambda p: np.clip(p + 0.18 * np.sin(3 * np.pi * p), 1e-6, 1 - 1e-6), "B: non-monotone distortion")
# A: sigmoidal distortion
#   n_cal      raw   sigmoid  isotonic   winner
#     200   0.2296    0.2009    0.2058   sigmoid
#    1000   0.2300    0.2001    0.2021   sigmoid
#    5000   0.2294    0.1996    0.2002   sigmoid
#   20000   0.2304    0.2001    0.2003   sigmoid
#
# B: non-monotone distortion
#   n_cal      raw   sigmoid  isotonic   winner
#     200   0.2159    0.2223    0.2131   isotonic
#    1000   0.2160    0.2223    0.2090   isotonic
#    5000   0.2159    0.2216    0.2075   isotonic
#   20000   0.2161    0.2219    0.2074   isotonic

2. Calibrating on the training set. Fit a calibrator on the model's own training predictions and on held-out data, and compare both under Brier and log loss.

You should get: worse test scores after fitting calibration on training predictions, with a larger relative increase in log loss.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.calibration import CalibratedClassifierCV
from sklearn.frozen import FrozenEstimator
from sklearn.metrics import brier_score_loss, log_loss
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=40000, n_features=12, n_informative=6,
                           class_sep=0.7, random_state=0)
Xtr, Xrest, ytr, yrest = train_test_split(X, y, test_size=0.6, random_state=0)
Xpool, Xte, ypool, yte = train_test_split(Xrest, yrest, test_size=0.5, random_state=0)
rf = RandomForestClassifier(n_estimators=300, random_state=0, n_jobs=-1).fit(Xtr, ytr)

for label, Xc, yc in (("training data", Xtr, ytr), ("held-out data", Xpool, ypool)):
    c = CalibratedClassifierCV(FrozenEstimator(rf), method="isotonic").fit(Xc, yc)
    p = c.predict_proba(Xte)[:, 1]
    print(f"calibrated on {label:14s} test Brier {brier_score_loss(yte, p):.4f}"
          f"  test logloss {log_loss(yte, p):.4f}")
# calibrated on training data  test Brier 0.0493  test logloss 1.4251
# calibrated on held-out data  test Brier 0.0429  test logloss 0.1551

3. Calibration without discrimination. Decompose the Brier score into reliability, resolution, and uncertainty, and compare a constant forecast, coarse probability estimates, and probabilities with a distorted scale.

You should get: a calibrated constant with no discrimination, and rescaled probabilities that rank well but have worse Brier than that constant. Distinguish original and binned forecasts.

import numpy as np
from sklearn.metrics import brier_score_loss, roc_auc_score

rng = np.random.default_rng(0)
n = 200_000
p_true = rng.beta(2, 5, n)
y = (rng.random(n) < p_true).astype(int)

def decompose(y, p, n_bins=20):
    edges = np.quantile(p, np.linspace(0, 1, n_bins + 1)); edges[-1] += 1e-9
    idx = np.clip(np.searchsorted(edges, p, side="right") - 1, 0, n_bins - 1)
    base = y.mean(); rel = res = 0.0
    grouped = np.empty_like(p)
    for b in range(n_bins):
        m = idx == b
        if not m.any(): continue
        w = m.mean()
        grouped[m] = p[m].mean()
        rel += w * (p[m].mean() - y[m].mean()) ** 2      # binned reliability
        res += w * (y[m].mean() - base) ** 2             # binned resolution
    unc = base * (1 - base)
    binned_brier = brier_score_loss(y, grouped)
    assert np.isclose(binned_brier, rel - res + unc)
    return rel, res, unc, binned_brier

models = {"true probabilities": p_true,
          "known base rate": np.full(n, 2 / 7),
          "coarse estimates": np.where(p_true > 0.28, 0.42, 0.20),
          "rescaled probabilities": np.clip(p_true * 2.4, 0, 1)}
print(f"{'model':28s} {'Brier':>7} {'binned':>7} {'AUC':>7} {'reliab':>8} {'resol':>8} {'uncert':>8}")
for name, p in models.items():
    rel, res, unc, binned = decompose(y, p)
    print(f"{name:28s} {brier_score_loss(y, p):7.4f} {binned:7.4f} {roc_auc_score(y, p):7.4f}"
          f" {rel:8.4f} {res:8.4f} {unc:8.4f}")
# model                          Brier  binned     AUC   reliab    resol   uncert
# true probabilities            0.1793  0.1795  0.7186   0.0000   0.0248   0.2044
# known base rate               0.2044  0.2044  0.5000   0.0000   0.0000   0.2044
# coarse estimates              0.1880  0.1880  0.6592   0.0007   0.0170   0.2044
# rescaled probabilities        0.3211  0.3210  0.7135   0.1402   0.0235   0.2044