{ "cells": [ { "cell_type": "markdown", "id": "7fb27b941602401d91542211134fc71a", "metadata": {}, "source": [ "# Evaluating a Calibrator Honestly\n", "\n", "Two questions this notebook answers:\n", "\n", "1. **How do you measure calibration without fooling yourself?** Scoring a calibrator on\n", " its own training data does not merely flatter it — for an isotonic-family calibrator\n", " it reports perfect calibration *by construction*.\n", "2. **When a model gets worse, which part broke?** A single score says \"worse\". The CORP\n", " decomposition says how much is miscalibration you can fix and how much is lost\n", " discrimination you cannot.\n", "\n", "Both rest on the same machinery: isotonic regression via the pool-adjacent-violators\n", "algorithm, which is what `calibre` is built on throughout.\n", "\n", "Reference: Dimitriadis, Gneiting & Jordan (2021), *Stable reliability diagrams for\n", "probabilistic classifiers*, PNAS 118(8)." ] }, { "cell_type": "code", "execution_count": null, "id": "acae54e37e7d407bbb7b55eff062a284", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from calibre import (\n", " IsotonicCalibrator,\n", " confidence_bands,\n", " consistency_bands,\n", " corp_reliability,\n", " cross_val_calibrate,\n", " debiased_calibration_error,\n", " score_decomposition,\n", " sweep_calibration_error,\n", ")\n", "from calibre.metrics import expected_calibration_error\n", "\n", "rng = np.random.default_rng(20260731)" ] }, { "cell_type": "markdown", "id": "9a63283cbaf04dbcab1f6479b197f3a8", "metadata": {}, "source": [ "## 1. A reliability diagram with no bins to choose\n", "\n", "The classic reliability diagram makes you pick a bin count, and the picture changes with\n", "the choice. Below, the same forecasts are binned three ways — the apparent calibration\n", "swings with the bin count, which is an artifact of the analyst's choice rather than a\n", "property of the model.\n", "\n", "CORP removes the choice: isotonic regression decides the number and position of the flat\n", "segments, optimally and automatically." ] }, { "cell_type": "code", "execution_count": null, "id": "8dd0d8092fe74a7c96281538738b07e2", "metadata": {}, "outputs": [], "source": [ "n = 1000\n", "truth = rng.uniform(0, 1, n)\n", "labels = rng.binomial(1, truth).astype(float)\n", "# A mildly overconfident forecaster.\n", "forecasts = np.clip(1.4 * (truth - 0.5) + 0.5, 0.001, 0.999)\n", "\n", "\n", "def binned_curve(x, y, n_bins):\n", " edges = np.quantile(x, np.linspace(0, 1, n_bins + 1))\n", " idx = np.clip(np.digitize(x, edges) - 1, 0, n_bins - 1)\n", " xs, ys = [], []\n", " for b in range(n_bins):\n", " m = idx == b\n", " if m.sum():\n", " xs.append(x[m].mean())\n", " ys.append(y[m].mean())\n", " return np.array(xs), np.array(ys)\n", "\n", "\n", "fig, axes = plt.subplots(1, 4, figsize=(16, 4), sharex=True, sharey=True)\n", "for ax, n_bins in zip(axes[:3], (5, 10, 20)):\n", " xs, ys = binned_curve(forecasts, labels, n_bins)\n", " ax.plot([0, 1], [0, 1], \"k--\", lw=1)\n", " ax.plot(xs, ys, \"o-\", color=\"crimson\")\n", " ax.set_title(f\"binned, {n_bins} bins\")\n", " ax.set_xlabel(\"forecast\")\n", "\n", "diagram = corp_reliability(forecasts, labels)\n", "axes[3].plot([0, 1], [0, 1], \"k--\", lw=1)\n", "axes[3].step(diagram.x, diagram.cep, color=\"crimson\", where=\"post\")\n", "axes[3].set_title(\"CORP (no bin count)\")\n", "axes[3].set_xlabel(\"forecast\")\n", "axes[0].set_ylabel(\"observed frequency\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "72eea5119410473aa328ad9291626812", "metadata": {}, "source": [ "## 2. Is the deviation real, or is it noise?\n", "\n", "A curve off the diagonal means nothing without knowing how far a *calibrated* forecaster\n", "would stray by chance. Consistency bands answer exactly that: outcomes are redrawn as\n", "`y* ~ Bernoulli(x)` — taking the forecasts at face value — and the diagram refit. The\n", "bands sit around the diagonal, and a curve leaving them is the analogue of a small\n", "p-value.\n", "\n", "Confidence bands answer the other question, clustering around the estimate instead." ] }, { "cell_type": "code", "execution_count": null, "id": "8edb47106e1a46a883d545849b8ab81b", "metadata": {}, "outputs": [], "source": [ "band = consistency_bands(forecasts, labels, level=0.9, n_resamples=400, random_state=0)\n", "conf = confidence_bands(forecasts, labels, level=0.9, n_resamples=400, random_state=0)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)\n", "for ax, (b, title) in zip(\n", " axes,\n", " (\n", " (band, \"consistency (around the diagonal)\"),\n", " (conf, \"confidence (around the estimate)\"),\n", " ),\n", "):\n", " ax.fill_between(b[\"x\"], b[\"lower\"], b[\"upper\"], alpha=0.3, color=\"steelblue\")\n", " ax.plot([0, 1], [0, 1], \"k--\", lw=1)\n", " ax.step(diagram.x, diagram.cep, color=\"crimson\", where=\"post\")\n", " ax.set_title(title)\n", " ax.set_xlabel(\"forecast\")\n", "axes[0].set_ylabel(\"conditional event probability\")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "outside = np.mean((diagram.cep < band[\"lower\"]) | (diagram.cep > band[\"upper\"]))\n", "print(f\"fraction of the curve outside the 90% consistency band: {outside:.1%}\")" ] }, { "cell_type": "markdown", "id": "10185d26023b46108eb7d9f57d49d2b3", "metadata": {}, "source": [ "## 3. Where did the score go?\n", "\n", "`score_decomposition` splits any proper score into three parts:\n", "\n", " mean_score = MCB - DSC + UNC\n", "\n", "- **MCB** (miscalibration) — what recalibration would save you. Fixable.\n", "- **DSC** (discrimination) — what your scores buy over always predicting the base rate.\n", "- **UNC** (uncertainty) — the difficulty of the problem. Nobody can change it.\n", "\n", "`MCB` and `DSC` are non-negative by construction, and the identity is exact." ] }, { "cell_type": "code", "execution_count": null, "id": "8763a12b2bbd4a93a75aff182afb95dc", "metadata": {}, "outputs": [], "source": [ "variants = {\n", " \"honest\": truth,\n", " \"overconfident\": np.clip(1.6 * (truth - 0.5) + 0.5, 0, 1),\n", " \"underconfident\": 0.5 * (truth - 0.5) + 0.5,\n", " \"noise-added\": np.clip(truth + rng.normal(0, 0.15, n), 0, 1),\n", "}\n", "\n", "print(f\"{'forecaster':16s} {'Brier':>8s} {'MCB':>8s} {'DSC':>8s} {'UNC':>8s}\")\n", "for name, x in variants.items():\n", " d = score_decomposition(x, labels)\n", " print(\n", " f\"{name:16s} {d['mean_score']:8.4f} {d['MCB']:8.4f} \"\n", " f\"{d['DSC']:8.4f} {d['UNC']:8.4f}\"\n", " )" ] }, { "cell_type": "markdown", "id": "7623eae2785240b9bd12b16a66d81610", "metadata": {}, "source": [ "Read the table by column, not by row.\n", "\n", "`UNC` is identical everywhere — it depends only on the outcomes, so it is the same\n", "problem in every row. The forecasters differ in the other two, and they differ in\n", "*different ways*: a monotone squeeze damages `MCB` while leaving `DSC` largely intact,\n", "because squeezing preserves the ranking. Adding noise damages `DSC`, because it destroys\n", "ranking information that no recalibration can recover.\n", "\n", "That distinction is the practical payoff. A high `MCB` is a call to recalibrate. A low\n", "`DSC` is a call to build a better model." ] }, { "cell_type": "markdown", "id": "7cdc8c89c7104fffa095e18ddfef8986", "metadata": {}, "source": [ "## 4. The trap: never evaluate on the training data\n", "\n", "For an isotonic-family calibrator this is not a matter of degree. The calibrator and the\n", "CORP diagnostic are the *same* PAV projection, and PAV is idempotent — so in-sample\n", "`MCB` is exactly zero however badly the model generalises. It cannot detect\n", "miscalibration even in principle.\n", "\n", "`cross_val_calibrate` returns out-of-fold probabilities: each one from a model that\n", "never saw that observation." ] }, { "cell_type": "code", "execution_count": null, "id": "b118ea5561624da68c537baed56e602f", "metadata": {}, "outputs": [], "source": [ "in_sample = IsotonicCalibrator().fit(forecasts, labels).transform(forecasts)\n", "out_of_fold = cross_val_calibrate(IsotonicCalibrator(), forecasts, labels, cv=5)\n", "\n", "for name, values in ((\"in-sample\", in_sample), (\"out-of-fold\", out_of_fold)):\n", " d = score_decomposition(values, labels)\n", " print(f\"{name:12s} MCB {d['MCB']:.6f} Brier {d['mean_score']:.4f}\")" ] }, { "cell_type": "markdown", "id": "938c804e27f84196a10c8828c723f798", "metadata": {}, "source": [ "The in-sample `MCB` is zero to machine precision. It would be zero for a calibrator\n", "fit to pure noise too. Any number you intend to believe should be computed out-of-fold." ] }, { "cell_type": "markdown", "id": "504fb2a444614c0babb325280ed9130a", "metadata": {}, "source": [ "## 5. Binned calibration error is biased\n", "\n", "If you do want a single calibration-error number, know that the plugin binned estimator\n", "reports error that is not there: part of every bin's gap is sampling noise in the label\n", "mean. The bias grows with the bin count — exactly when you wanted a finer picture." ] }, { "cell_type": "code", "execution_count": null, "id": "59bbdb311c014d738909a11f9e486628", "metadata": {}, "outputs": [], "source": [ "calibrated = rng.uniform(0, 1, 4000)\n", "calibrated_y = rng.binomial(1, calibrated).astype(float) # true error is 0\n", "\n", "rows = [\n", " (\n", " n_bins,\n", " expected_calibration_error(calibrated_y, calibrated, n_bins=n_bins),\n", " debiased_calibration_error(calibrated_y, calibrated, n_bins=n_bins),\n", " )\n", " for n_bins in (5, 10, 20, 50)\n", "]\n", "\n", "print(f\"{'bins':>5s} {'plugin':>9s} {'debiased':>9s} (true error is 0)\")\n", "for n_bins, plugin, deb in rows:\n", " print(f\"{n_bins:5d} {plugin:9.4f} {deb:9.4f}\")\n", "\n", "print()\n", "print(\n", " f\"sweep (chooses its own bin count): \"\n", " f\"{sweep_calibration_error(calibrated_y, calibrated):.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "b43b363d81ae4b689946ece5c682cd59", "metadata": {}, "source": [ "The plugin estimate climbs with the bin count on data with no miscalibration at all.\n", "The debiased estimator stays near zero.\n", "\n", "`sweep_calibration_error` attacks the same problem from the other side: it adds bins\n", "while the calibration curve stays monotone and stops when it doesn't, on the reasoning\n", "that non-monotonicity is the signal the bins have become fine enough to read noise.\n", "\n", "Both use equal-mass bins, and neither ever splits a group of tied predictions across a\n", "bin boundary — which matters more than it sounds, since clipping a forecast into [0, 1]\n", "routinely puts hundreds of observations on a single value." ] }, { "cell_type": "markdown", "id": "8a65eabff63a45729fe45fb5ade58bdc", "metadata": {}, "source": [ "## Summary\n", "\n", "- Use `corp_reliability` for a reliability diagram with nothing to tune, and\n", " `consistency_bands` to see whether the deviation is real.\n", "- Use `score_decomposition` to separate the miscalibration you can fix from the\n", " discrimination you cannot.\n", "- Compute all of it on `cross_val_calibrate` output. In-sample calibration error for an\n", " isotonic-family calibrator is identically zero and tells you nothing.\n", "- If you need a scalar, prefer `debiased_calibration_error` or\n", " `sweep_calibration_error` to the plugin ECE.\n", "\n", "These numbers are pinned against R's `reliabilitydiag` to machine precision in\n", "`tests/test_r_reference.py`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.12" } }, "nbformat": 4, "nbformat_minor": 5 }