Probability Calibration › Implementation in NumPy, SciPy, and Scikit-learn
Machine Learning from ScratchAll chapters Probability Calibration › Implementation in NumPy, SciPy, and Scikit-learnSection view · Sections in this chapter6 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
We now build the calibration workflow from primitives and compare it with standard library implementations. The implementation follows the logic of the formulas: fit the base classifier on training data, fit each calibrator on a separate calibration partition, and evaluate the results on an untouched test partition. We fit a Gaussian naive Bayes classifier. To make a controlled scale error visible, we then multiply its predicted log odds by two and convert them back to probabilities. This deliberately sharpens the scores without adding feature information; it is a synthetic demonstration, not a claim that this classifier is always overconfident.
The next cells import the libraries and set the data-generation seed, then generate and split the data. The calibration and test predictions use rows held out from base-model fitting.
PythonIn [7]Expand
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from book_plot_style import (
PALETTE,
polish_axes,
use_book_style,
)
from scipy.optimize import minimize
from scipy.special import expit
from sklearn.datasets import make_classification
from sklearn.isotonic import IsotonicRegression
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
RANDOM_STATE = 42
EPS = 1e-12
SCORE_SHARPNESS = 2.0PythonIn [8]Expand
X, y = make_classification(
n_samples=12000,
n_features=12,
n_informative=6,
n_redundant=2,
n_clusters_per_class=2,
class_sep=0.9,
flip_y=0.04,
random_state=RANDOM_STATE,
)
X_train, X_temp, y_train, y_temp = train_test_split(
X, y, test_size=0.60, stratify=y, random_state=RANDOM_STATE
)
X_cal, X_test, y_cal, y_test = train_test_split(
X_temp, y_temp, test_size=0.50, stratify=y_temp, random_state=RANDOM_STATE
)
base = GaussianNB()
base.fit(X_train, y_train)
p_cal_base = base.predict_proba(X_cal)[:, 1]
p_test_base = base.predict_proba(X_test)[:, 1]Shared input checks reject invalid probabilities, nonbinary outcomes and mismatched row counts before arithmetic can silently broadcast them. We then define numerical clipping, probability scores and reliability-bin summaries. The metric table is then displayed from variables so the comparison is computed rather than typed.
PythonIn [9]Expand
def probability_vector(p):
p = np.asarray(p, dtype=float)
if (
p.ndim != 1
or p.size == 0
or not np.isfinite(p).all()
or np.any((p < 0) | (p > 1))
):
raise ValueError(
"Probabilities must be a nonempty finite vector in [0, 1]."
)
return p
def binary_pairs(y, p):
p = probability_vector(p)
y = np.asarray(y, dtype=float)
if (
y.ndim != 1
or y.shape != p.shape
or not np.isfinite(y).all()
or not np.isin(y, [0, 1]).all()
):
raise ValueError(
"Binary outcomes must have the same one-dimensional shape as probabilities."
)
return y, p
def clip_prob(p, eps=EPS):
p = probability_vector(p)
if not np.isfinite(eps) or not 0 < eps < 0.5:
raise ValueError("eps must be finite and between zero and one half.")
return np.clip(p, eps, 1.0 - eps)
def brier_score(y, p):
y, p = binary_pairs(y, p)
return float(np.mean((p - y) ** 2))
def log_loss_score(y, p):
y, p = binary_pairs(y, p)
p = clip_prob(p)
return float(-np.mean(y * np.log(p) + (1.0 - y) * np.log1p(-p)))
def reliability_table(y, p, n_bins=10, strategy="quantile"):
y, p = binary_pairs(y, p)
if (
not isinstance(n_bins, (int, np.integer))
or isinstance(n_bins, bool)
or n_bins < 1
):
raise ValueError("n_bins must be a positive integer.")
if strategy == "quantile":
edges = np.quantile(p, np.linspace(0.0, 1.0, n_bins + 1))
edges = np.unique(edges)
elif strategy == "uniform":
edges = np.linspace(0.0, 1.0, n_bins + 1)
else:
raise ValueError("strategy must be 'quantile' or 'uniform'")
if len(edges) < 2:
edges = np.array([0.0, 1.0])
bin_idx = np.digitize(p, edges[1:-1], right=False)
rows = []
for b in range(len(edges) - 1):
mask = bin_idx == b
if not np.any(mask):
continue
rows.append(
{
"bin": b,
"n": int(mask.sum()),
"mean_pred": float(p[mask].mean()),
"observed": float(y[mask].mean()),
}
)
return pd.DataFrame(rows)
def expected_calibration_error(y, p, n_bins=10, strategy="quantile"):
tab = reliability_table(y, p, n_bins=n_bins, strategy=strategy)
n = len(y)
return float(
np.sum((tab["n"] / n) * np.abs(tab["observed"] - tab["mean_pred"]))
)PythonIn [10]Expand
def logit_clipped(p, eps=EPS):
p = clip_prob(p, eps)
return np.log(p / (1.0 - p))
def fit_sigmoid(p, y):
y, p = binary_pairs(y, p)
if np.unique(y).size != 2:
raise ValueError("Sigmoid fitting requires both outcome classes.")
z = logit_clipped(p)
def objective(theta):
a, b = theta
eta = a * z + b
return np.mean(np.logaddexp(0.0, eta) - y * eta)
def gradient(theta):
a, b = theta
error = expit(a * z + b) - y
return np.array([np.mean(error * z), np.mean(error)])
result = minimize(
objective,
x0=np.array([1.0, 0.0]),
jac=gradient,
method="BFGS",
options={"gtol": 1e-7, "maxiter": 1000},
)
if (
not result.success
or not np.isfinite(result.x).all()
or not np.isfinite(result.fun)
):
raise RuntimeError(
"Sigmoid optimization did not return a finite successful fit."
)
return result.x
def predict_sigmoid(p, theta):
theta = np.asarray(theta, dtype=float)
if theta.shape != (2,) or not np.isfinite(theta).all():
raise ValueError("theta must contain two finite coefficients.")
a, b = theta
return expit(a * logit_clipped(p) + b)The isotonic calibrator uses the pool adjacent violators algorithm. We first aggregate tied scores, then merge adjacent blocks whenever their observed rates decrease. The fitted block means are the calibrated probabilities. Prediction uses linear interpolation between the unique score thresholds and holds the boundary values constant outside the observed range.
PythonIn [11]Expand
def fit_isotonic(p, y):
y, p = binary_pairs(y, p)
order = np.argsort(p)
p_sorted = p[order]
y_sorted = y[order]
unique_p, inv = np.unique(p_sorted, return_inverse=True)
sum_y = np.bincount(inv, weights=y_sorted)
counts = np.bincount(inv)
blocks = []
for i in range(len(unique_p)):
blocks.append([sum_y[i], counts[i], i, i])
while len(blocks) >= 2:
prev = blocks[-2]
curr = blocks[-1]
if prev[0] / prev[1] <= curr[0] / curr[1]:
break
blocks.pop()
blocks.pop()
merged = [prev[0] + curr[0], prev[1] + curr[1], prev[2], curr[3]]
blocks.append(merged)
fitted = np.empty(len(unique_p), dtype=float)
for s, w, start, end in blocks:
fitted[start : end + 1] = s / w
return unique_p, fitted
def predict_isotonic(p, iso):
p = probability_vector(p)
unique_p, fitted = iso
return np.interp(p, unique_p, fitted, left=fitted[0], right=fitted[-1])We sharpen the base scores on the logit scale, then fit both calibrators using only the calibration labels and evaluate on the test partition. The Brier score and log loss are computed from the test labels only. AUC is included to compare ranking with probability quality.
PythonIn [12]Expand
p_cal_raw = expit(SCORE_SHARPNESS * logit_clipped(p_cal_base))
p_test_raw = expit(SCORE_SHARPNESS * logit_clipped(p_test_base))
theta = fit_sigmoid(p_cal_raw, y_cal)
iso = fit_isotonic(p_cal_raw, y_cal)
p_test_sig = predict_sigmoid(p_test_raw, theta)
p_test_iso = predict_isotonic(p_test_raw, iso)
metrics = pd.DataFrame(
{
"model": ["Raw", "Sigmoid", "Isotonic"],
"Brier": [
brier_score(y_test, p_test_raw),
brier_score(y_test, p_test_sig),
brier_score(y_test, p_test_iso),
],
"Log loss": [
log_loss_score(y_test, p_test_raw),
log_loss_score(y_test, p_test_sig),
log_loss_score(y_test, p_test_iso),
],
"AUC": [
roc_auc_score(y_test, p_test_raw),
roc_auc_score(y_test, p_test_sig),
roc_auc_score(y_test, p_test_iso),
],
}
)PythonIn [13]Expand
def summarize_metrics(table):
return table.round(4)
def format_calibration_summary(theta, y_cal, y_test):
return (
f"Sigmoid calibration: a = {theta[0]:.3f}, b = {theta[1]:.3f}\n"
f"Calibration set size: {len(y_cal)}; test set size: {len(y_test)}"
)| model | Brier | Log loss | AUC | |
|---|---|---|---|---|
| 0 | Raw | 0.1609 | 0.5497 | 0.8691 |
| 1 | Sigmoid | 0.1473 | 0.4548 | 0.8691 |
| 2 | Isotonic | 0.1478 | 0.4560 | 0.8670 |
The following display reports the fitted sigmoid coefficients and the actual partition sizes.
The fitted sigmoid and isotonic maps estimate a correction to the probability scale. The next figure shows how each map transforms raw predicted probabilities on the test range, with the identity line as a reference for no change.
PythonIn [15]Expand
p_grid = np.linspace(float(p_test_raw.min()), float(p_test_raw.max()), 300)
sig_grid = predict_sigmoid(p_grid, theta)
iso_grid = predict_isotonic(p_grid, iso)VisualizationOut [16]Expand
ConsoleOut [17]Expand
Sigmoid calibration: a = 0.483, b = -0.117 Calibration set size: 3600; test set size: 3600
The fitted sigmoid coefficients describe the change in log odds. With zero intercept, a positive slope below one moves probabilities toward one half, while a slope above one moves them toward the endpoints. A nonzero intercept also shifts the curve, so the slope alone does not describe the adjustment. The identity mapping has slope one and intercept zero.
The binning sensitivity table repeats the expected calibration error for several bin counts. The numbers are not expected to be identical. A diagnostic that changes with the binning scheme should not be treated as a definitive model score.
PythonIn [18]Expand
ece_rows = []
for name, p in [
("Raw", p_test_raw),
("Sigmoid", p_test_sig),
("Isotonic", p_test_iso),
]:
for bins in [5, 10, 20]:
ece_rows.append(
{
"model": name,
"bins": bins,
"ECE": expected_calibration_error(
y_test, p, n_bins=bins, strategy="quantile"
),
}
)
ece_table = pd.DataFrame(ece_rows)| bins | 5 | 10 | 20 |
|---|---|---|---|
| model | |||
| Isotonic | 0.0067 | 0.0151 | 0.0202 |
| Raw | 0.1049 | 0.1049 | 0.1049 |
| Sigmoid | 0.0193 | 0.0201 | 0.0257 |
We can compare our implementation with scikit-learn. The reference sigmoid fit is a logistic regression on the clipped logit, and the reference isotonic fit uses IsotonicRegression with boundary clipping. Small differences can arise from numerical optimization and interpolation, but the calibrated probabilities should be close. The maximum absolute differences check agreement on these test scores. Optimization tolerances and the reference logistic regression penalty can leave a small gap; agreement on this dataset is a useful implementation check.
PythonIn [20]Expand
iso_ref = IsotonicRegression(out_of_bounds="clip").fit(p_cal_raw, y_cal)
p_test_iso_ref = iso_ref.predict(p_test_raw)
z_cal = logit_clipped(p_cal_raw).reshape(-1, 1)
z_test = logit_clipped(p_test_raw).reshape(-1, 1)
lr = LogisticRegression(C=1e6, solver="lbfgs", max_iter=1000)
lr.fit(z_cal, y_cal)
p_test_sig_ref = lr.predict_proba(z_test)[:, 1]ConsoleOut [21]Expand
Max absolute difference, sigmoid: 0.000044 Max absolute difference, isotonic: 0.000000
Key Parameters
Below are the main parameters that affect how the calibrators work and how reliable the resulting probabilities are.
Calibration-set size and support across the score range are important practical constraints. Sigmoid recalibration has only two parameters and can work with a relatively modest calibration sample, but isotonic recalibration can fit many blocks and needs more data to avoid overfitting. There is no universal threshold that works for every problem; the relevant question is whether the calibration set contains enough positives and negatives in the score regions that matter for decisions.
The choice is a flexibility tradeoff. Sigmoid fits a smooth two-parameter curve, which can miss a more irregular calibration relationship. Isotonic fits a nondecreasing relationship with plateaus that can create ties; it can overfit a small calibration sample. Check support across the score range and compare methods using validation data that remains separate from final evaluation.
The reliability-diagram parameters are the number of bins and the binning strategy. Quantile bins put roughly the same number of observations in each bin, which stabilizes the observed rates but can widen bins in sparse probability regions. Uniform bins keep the probability axis evenly spaced but can produce noisy extreme bins. The ECE depends on both choices, so we report it alongside Brier score and log loss rather than in isolation.
EPS clips only logits and log-loss evaluation, so the computed log loss is that of the clipped forecast. The ideal log loss can be infinite at a confidently wrong endpoint. Brier scoring uses the original validated probabilities. SCORE_SHARPNESS is fixed at two for this controlled example; it is not selected using the test labels.
Key Methods
The following are the main functions we use for calibrating and evaluating probabilities.
-
fit_sigmoidandpredict_sigmoidimplement Platt scaling on the logit scale. -
fit_isotonicandpredict_isotonicimplement the pool adjacent violators algorithm and monotone interpolation. -
brier_scorecomputes the exact binary squared-error score, whilelog_loss_scoreevaluates the numerically clipped forecast on the test set. -
reliability_tablebins predictions for reliability diagrams. -
expected_calibration_errorsummarizes the binned gaps. -
roc_auc_scoremeasures discrimination and is used to verify that calibration has not been mistaken for a ranking improvement.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!