PD calibration walkthrough¶
A full probcal cycle on a synthetic credit-risk portfolio: diagnose → select → fit → re-anchor → backtest → translate thresholds. Every step links to the concept chapter that explains it.
The portfolio is deliberately miscalibrated the way low-default portfolios are in practice: a spread error plus an asymmetric low-PD tail distortion, at a ~3% event rate.
In [1]:
Copied!
import numpy as np
from probcal import (
BetaCalibrator,
CalibratorSelector,
LogitOffset,
make_pd_portfolio,
)
from probcal.metrics import calibration_guardrails, evaluate, jeffreys_grade_test
port = make_pd_portfolio(n=4000, random_state=42)
s, y = port.scores, port.y
print(f'portfolio: n={len(s)}, defaults={int(y.sum())}, mean score={s.mean():.4f}, '
f'realized rate={y.mean():.4f}')
import numpy as np
from probcal import (
BetaCalibrator,
CalibratorSelector,
LogitOffset,
make_pd_portfolio,
)
from probcal.metrics import calibration_guardrails, evaluate, jeffreys_grade_test
port = make_pd_portfolio(n=4000, random_state=42)
s, y = port.scores, port.y
print(f'portfolio: n={len(s)}, defaults={int(y.sum())}, mean score={s.mean():.4f}, '
f'realized rate={y.mean():.4f}')
portfolio: n=4000, defaults=114, mean score=0.0613, realized rate=0.0285
1. Diagnose¶
The guardrail triplet (slope, intercept, Spiegelhalter) localizes the defect; see Metrics and tests.
In [2]:
Copied!
g = calibration_guardrails(y, s)
print(f'slope {g.slope:+.3f} ok={g.slope_ok}')
print(f'intercept {g.intercept:+.3f} ok={g.intercept_ok}')
print(f'spiegelhalter p {g.spiegelhalter_p:.4f} ok={g.spiegelhalter_ok}')
g = calibration_guardrails(y, s)
print(f'slope {g.slope:+.3f} ok={g.slope_ok}')
print(f'intercept {g.intercept:+.3f} ok={g.intercept_ok}')
print(f'spiegelhalter p {g.spiegelhalter_p:.4f} ok={g.spiegelhalter_ok}')
slope +0.725 ok=False intercept -0.856 ok=False spiegelhalter p 0.0000 ok=False
The logit-scale reliability diagram is the recommended view for low-PD work — the low-probability region gets the resolution the decisions need (Visualization).
In [3]:
Copied!
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from probcal import reliability_binned, reliability_loess
from probcal.plots import plot_reliability
ax = plot_reliability(reliability_binned(y, s), smooth=reliability_loess(y, s), scale='logit')
plt.show()
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from probcal import reliability_binned, reliability_loess
from probcal.plots import plot_reliability
ax = plot_reliability(reliability_binned(y, s), smooth=reliability_loess(y, s), scale='logit')
plt.show()
/tmp/ipykernel_652022/743357854.py:9: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown plt.show()
2. Select¶
CalibratorSelector compares candidates by out-of-fold log loss — never on their own fitting data (Automatic selection). A restricted menu keeps this demo fast.
In [4]:
Copied!
from probcal import CenteredIsotonicCalibrator, PlattCalibrator, TemperatureCalibrator
sel = CalibratorSelector(
candidates={
'platt': PlattCalibrator(),
'temperature': TemperatureCalibrator(),
'beta_abm': BetaCalibrator(),
'cir': CenteredIsotonicCalibrator(),
},
cv=5,
random_state=0,
).fit(s, y)
print(sel.report_)
from probcal import CenteredIsotonicCalibrator, PlattCalibrator, TemperatureCalibrator
sel = CalibratorSelector(
candidates={
'platt': PlattCalibrator(),
'temperature': TemperatureCalibrator(),
'beta_abm': BetaCalibrator(),
'cir': CenteredIsotonicCalibrator(),
},
cv=5,
random_state=0,
).fit(s, y)
print(sel.report_)
SelectionReport (criterion: log_loss) method log_loss sd guardrails chosen ----------- -------- ---------- ---------- ------ beta_abm 0.120866 0.00189116 True platt 0.12102 0.00175219 True * temperature 0.126539 0.00343572 False cir 0.126635 0.0130849 False
3. Fit and interpret¶
The winner is refitted on the full calibration set; interpret() is the auditable reading of every fitted parameter (Parametric methods).
In [5]:
Copied!
cal = sel.best_calibrator_
print(sel.best_name_)
print()
print(cal.interpret())
cal = sel.best_calibrator_
print(sel.best_name_)
print()
print(cal.interpret())
platt Interpretation[PlattCalibrator] parameter value --------- -------- a 0.718393 b -1.50335 - slope a = 0.718 < 1: scores were overconfident (too spread out); predictions are shrunk toward the base rate - intercept b = -1.503: base-rate (calibration-in-the-large) shift of -1.503 log-odds, odds factor 0.222 - identity map corresponds to (a, b) = (1, 0)
In [6]:
Copied!
p_cal = cal.predict_proba(s)
rep = evaluate(y, p_cal, n_boot=60, seed=0)
print(rep)
p_cal = cal.predict_proba(s)
rep = evaluate(y, p_cal, n_boot=60, seed=0)
print(rep)
MetricReport metric value ci_low ci_high --------------- ----------- ----------- ---------- log_loss 0.120962 0.106428 0.134292 brier 0.0271773 0.0230925 0.0308998 brier_skill 0.0184353 0.0075793 0.0290465 ece 0.00718561 0.00722274 0.0158614 ece_debiased 0.00216233 0.00346722 0.0125181 mce 0.0142435 0.0171775 0.0451844 ece_sweep 0.00295916 0.00168338 0.00818723 smooth_ece 0.0159378 0.0180679 0.0222943 ecce_max 0.00181138 0.0020032 0.00601789 ecce_mean 0.000690679 0.000572931 0.00217352 ici 0.0040879 0.00296267 0.00890415 e50 0.00226302 0.00116732 0.00602503 e90 0.0127443 0.00478659 0.0256834 emax 0.127652 0.0449838 0.174722 spiegelhalter_z 0.0328449 -1.74397 1.55859 spiegelhalter_p 0.973798 0.0730885 0.963855 intercept -0.00015432 -0.18664 0.137798 slope 1.00957 0.775489 1.22777
In [7]:
Copied!
off = LogitOffset(target_mean=0.031).fit(p_cal)
p_final = off.transform(p_cal)
print(off.audit_report(y, p_cal))
off = LogitOffset(target_mean=0.031).fit(p_cal)
p_final = off.transform(p_cal)
print(off.audit_report(y, p_cal))
AuditReport(delta=+0.0885, odds factor 1.0925, fitted 2026-07-23T06:03:32+00:00) portfolio mean: 0.02850 -> 0.03100 slope: +1.010 -> +1.010 intercept: -0.000 -> -0.089 spiegelhalter p 0.974 -> 0.377 guardrails ok: True -> True
5. Per-grade backtest¶
Grades by calibrated-PD quantiles; the Jeffreys test is the ECB IRB formulation — one-sided and conservative (Metrics and tests).
In [8]:
Copied!
edges = np.quantile(p_final, [0.25, 0.5, 0.75])
grades = np.array(['G1', 'G2', 'G3', 'G4'])[np.searchsorted(edges, p_final)]
res = jeffreys_grade_test(y, p_final, grades)
for gname, n_g, k_g, pd_g, pv, light in zip(
res.grades, res.n, res.k, res.pd, res.p_value, res.light, strict=True
):
print(f'{gname}: n={n_g:5d} defaults={k_g:3d} PD={pd_g:.4f} '
f'p={pv:.3f} {light}')
edges = np.quantile(p_final, [0.25, 0.5, 0.75])
grades = np.array(['G1', 'G2', 'G3', 'G4'])[np.searchsorted(edges, p_final)]
res = jeffreys_grade_test(y, p_final, grades)
for gname, n_g, k_g, pd_g, pv, light in zip(
res.grades, res.n, res.k, res.pd, res.p_value, res.light, strict=True
):
print(f'{gname}: n={n_g:5d} defaults={k_g:3d} PD={pd_g:.4f} '
f'p={pv:.3f} {light}')
G1: n= 1000 defaults= 9 PD=0.0092 p=0.508 green G2: n= 1000 defaults= 15 PD=0.0186 p=0.800 green G3: n= 1000 defaults= 25 PD=0.0304 p=0.841 green G4: n= 1000 defaults= 65 PD=0.0657 p=0.529 green
6. Translate policy cutoffs to raw scores¶
Decisions are made on calibrated PD; deployed systems cut on raw scores. The generalized inverse translates policy into score space and refuses unattainable targets (Inverse maps).
In [9]:
Copied!
from probcal import calibrated_bands_to_raw
lo_s, hi_s = cal.interval_inverse(0.0, 0.02)
print(f"'approve below 2% PD' -> raw score <= {hi_s:.5f}")
masterscale = {'A': (0.0, 0.01), 'B': (0.01, 0.03), 'C': (0.03, 0.10), 'D': (0.10, 1.0)}
for grade, (lo_r, hi_r) in calibrated_bands_to_raw(cal, masterscale).items():
print(f'{grade}: raw scores [{lo_r:.5f}, {hi_r:.5f}]')
from probcal import calibrated_bands_to_raw
lo_s, hi_s = cal.interval_inverse(0.0, 0.02)
print(f"'approve below 2% PD' -> raw score <= {hi_s:.5f}")
masterscale = {'A': (0.0, 0.01), 'B': (0.01, 0.03), 'C': (0.03, 0.10), 'D': (0.10, 1.0)}
for grade, (lo_r, hi_r) in calibrated_bands_to_raw(cal, masterscale).items():
print(f'{grade}: raw scores [{lo_r:.5f}, {hi_r:.5f}]')
'approve below 2% PD' -> raw score <= 0.03473 A: raw scores [0.00000, 0.01334] B: raw scores [0.01334, 0.06031] C: raw scores [0.06031, 0.27570] D: raw scores [0.27570, 1.00000]
For a counterfactual engine the same inverse defines the raw target (interop recipe, shown not run — see the FAQ for the Target.probability trap):
lo_z, hi_z = cal.interval_inverse(0.0, 0.02, space='logit')
target = treecf.Target.raw(range=(lo_z, hi_z))
Pass buffer_logit=m to keep counterfactuals valid under future re-anchorings of magnitude up to m.