{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":13451,"datasetId":654585,"databundleVersionId":1188070},{"sourceType":"kernelVersion","sourceId":298700524},{"sourceType":"kernelVersion","sourceId":298979112}],"dockerImageVersionId":31287,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"c6ef1147","cell_type":"markdown","source":"# Notebook 6 — Confidence Calibration\n**RSNA Intracranial Hemorrhage Detection**\n\nA well-calibrated model outputs probabilities that match empirical frequencies.\nAn overconfident model may output p=0.95 for cases that are only 60% correct.\n\nThis notebook:\n1. Measures raw calibration (ECE, reliability diagram)\n2. Applies **Temperature Scaling** (parametric, fast)\n3. Applies **Isotonic Regression** (non-parametric, flexible)\n4. Compares all three and selects the best method\n5. Checks overconfidence (extreme probability shrinkage after scaling)\n6. Re-optimises the decision threshold after calibration\n7. Defines and validates the **3-band triage confidence system**\n8. Reports ECE at high-confidence band (prob > 0.8)\n9. Saves calibration parameters for use in the report generator\n\n### Calibration circularity caveat\n> **Note:** Temperature scaling and isotonic regression are both fitted and evaluated\n> on the same validation set. This means **reported post-calibration ECE is slightly\n> optimistic**. In a production setting, a held-out calibration set should be used.\n> However, with 100K+ validation samples, the optimism is minimal.\n>\n> **Isotonic regression is especially affected:** being non-parametric, it can\n> perfectly memorise the calibration mapping, yielding ECE ≈ 0 on the training set.\n> To get an honest estimate, we report **5-fold cross-validated isotonic ECE**.\n\n### Slice-level calibration caveat\n> **Note:** Calibration is performed at the **slice level**. In deployment, predictions\n> are typically aggregated to the **patient level** (e.g., via max-pooling). Patient-level\n> aggregation changes the probability distribution — max-pooling inflates probabilities,\n> so slice-level calibration does **not** directly transfer. A patient-level calibration\n> step would be needed in production.\n\n### Required inputs\n- NB02 output: `manifest.csv` + `cache/` NPY arrays + `normalization_stats.json`\n- NB03 output: `best_model.pth`, `checkpoint.pth`","metadata":{}},{"id":"9705182b","cell_type":"code","source":"# ── Config ────────────────────────────────────────────────────────────────\nimport os, json\nimport numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nimport torchvision.models as models\nimport torchvision.transforms as T\nfrom torch.utils.data import Dataset, DataLoader\nfrom sklearn.isotonic import IsotonicRegression\nfrom sklearn.metrics import roc_auc_score, brier_score_loss\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nimport seaborn as sns\nfrom tqdm.auto import tqdm\n\nDEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'\n\nCACHE_INPUT_DIR = '/kaggle/input/notebooks/harshitghosh/nb02eda'\nNPY_CACHE_DIR   = f'{CACHE_INPUT_DIR}/cache'\nMANIFEST_PATH   = f'{CACHE_INPUT_DIR}/manifest.csv'\nMODEL_PATH      = '/kaggle/input/notebooks/harshitghosh/03nbeda/best_model.pth'\nCHECKPOINT      = '/kaggle/input/notebooks/harshitghosh/03nbeda/checkpoint.pth'\n\nARCH        = 'efficientnet_b0'\nIMG_SIZE    = 256\nBATCH_SIZE  = 64\nNUM_WORKERS = 4\nSEED        = 42\nECE_BINS    = 15    # number of equal-width bins for ECE\n\n# ─── Load normalization stats ────────────────────────────────────────────\n_norm_path = os.path.join(CACHE_INPUT_DIR, 'normalization_stats.json')\nif os.path.exists(_norm_path):\n    with open(_norm_path) as f:\n        _norm = json.load(f)\n    MEAN = _norm['mean']\n    STD  = _norm['std']\n    print(f'Dataset normalization: mean={MEAN}, std={STD}')\nelse:\n    MEAN = [0.485, 0.456, 0.406]\n    STD  = [0.229, 0.224, 0.225]\n    print(f'Using ImageNet defaults: mean={MEAN}, std={STD}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:46:58.767836Z","iopub.execute_input":"2026-02-20T18:46:58.768571Z","iopub.status.idle":"2026-02-20T18:47:11.081253Z","shell.execute_reply.started":"2026-02-20T18:46:58.76854Z","shell.execute_reply":"2026-02-20T18:47:11.08051Z"}},"outputs":[],"execution_count":null},{"id":"8eae7fc5","cell_type":"code","source":"# ── 1. Load model & get raw logits on validation set ─────────────────────\ndef build_model(arch):\n    if arch == 'efficientnet_b0':\n        m = models.efficientnet_b0(weights=None)\n        m.classifier = nn.Sequential(nn.Dropout(0.3), nn.Linear(m.classifier[1].in_features, 1))\n    elif arch == 'resnet50':\n        m = models.resnet50(weights=None)\n        m.fc = nn.Sequential(nn.Dropout(0.3), nn.Linear(m.fc.in_features, 1))\n    return m\n\n\nclass ICHDataset(Dataset):\n    def __init__(self, df, npy_root, transform):\n        self.df = df.reset_index(drop=True)\n        self.npy_root = npy_root\n        self.transform = transform\n    def __len__(self): return len(self.df)\n    def __getitem__(self, idx):\n        row  = self.df.iloc[idx]\n        path = os.path.join(self.npy_root, f'{row[\"image_id\"]}.npy')\n        try:\n            img = np.load(path)  # uint8  [0, 255]\n        except Exception:\n            img = np.zeros((IMG_SIZE, IMG_SIZE, 3), np.uint8)\n        return self.transform(img), torch.tensor(float(row['any']), dtype=torch.float32)\n\n\nmodel = build_model(ARCH)\nmodel.load_state_dict(torch.load(MODEL_PATH, map_location=DEVICE))\nmodel = model.to(DEVICE).eval()\n\nckpt = torch.load(CHECKPOINT, map_location='cpu')\nBASE_THRESHOLD = ckpt.get('best_thresh', 0.5)\n\nmanifest = pd.read_csv(MANIFEST_PATH)\nval_df   = manifest[manifest['split'] == 'val'].reset_index(drop=True)\n\nval_transform = T.Compose([T.ToPILImage(), T.ToTensor(), T.Normalize(mean=MEAN, std=STD)])\nval_ds = ICHDataset(val_df, NPY_CACHE_DIR, val_transform)\nval_loader = DataLoader(val_ds, BATCH_SIZE, shuffle=False,\n                        num_workers=NUM_WORKERS, pin_memory=True)\n\n@torch.no_grad()\ndef collect_logits(model, loader):\n    logits_, labels_ = [], []\n    for imgs, lbls in tqdm(loader, desc='Collecting logits'):\n        with torch.amp.autocast(device_type='cuda'):\n            out = model(imgs.to(DEVICE)).squeeze(1).cpu().float()\n        logits_.append(out); labels_.append(lbls)\n    return torch.cat(logits_).numpy(), torch.cat(labels_).numpy()\n\n\nraw_logits, labels = collect_logits(model, val_loader)\nraw_probs = torch.sigmoid(torch.tensor(raw_logits)).numpy()\n\nprint(f'Collected logits for {len(labels):,} samples')\nprint(f'Raw AUC: {roc_auc_score(labels, raw_probs):.5f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:47:11.082553Z","iopub.execute_input":"2026-02-20T18:47:11.082954Z","iopub.status.idle":"2026-02-20T18:47:54.055771Z","shell.execute_reply.started":"2026-02-20T18:47:11.082929Z","shell.execute_reply":"2026-02-20T18:47:54.054846Z"}},"outputs":[],"execution_count":null},{"id":"bc0e1546","cell_type":"code","source":"# ── 2. ECE calculation ────────────────────────────────────────────────────\ndef expected_calibration_error(probs: np.ndarray,\n                                labels: np.ndarray,\n                                n_bins: int = 15) -> tuple:\n    \"\"\"\n    Compute ECE and return (ece, bin_confs, bin_accs, bin_sizes) for plotting.\n    \"\"\"\n    bins      = np.linspace(0.0, 1.0, n_bins + 1)\n    bin_ids   = np.digitize(probs, bins, right=True) - 1\n    bin_ids   = np.clip(bin_ids, 0, n_bins - 1)\n\n    bin_confs = np.zeros(n_bins)\n    bin_accs  = np.zeros(n_bins)\n    bin_sizes = np.zeros(n_bins)\n\n    for b in range(n_bins):\n        mask = bin_ids == b\n        if mask.sum() == 0:\n            continue\n        bin_sizes[b] = mask.sum()\n        bin_confs[b] = probs[mask].mean()\n        bin_accs[b]  = labels[mask].mean()\n\n    ece = (np.abs(bin_accs - bin_confs) * bin_sizes / len(probs)).sum()\n    return float(ece), bin_confs, bin_accs, bin_sizes\n\n\nraw_ece, raw_bc, raw_ba, raw_bs = expected_calibration_error(raw_probs, labels)\nprint(f'Raw ECE: {raw_ece:.5f}  (lower is better; 0 = perfect calibration)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:47:54.057544Z","iopub.execute_input":"2026-02-20T18:47:54.057856Z","iopub.status.idle":"2026-02-20T18:47:54.068221Z","shell.execute_reply.started":"2026-02-20T18:47:54.057822Z","shell.execute_reply":"2026-02-20T18:47:54.067641Z"}},"outputs":[],"execution_count":null},{"id":"c6555aba","cell_type":"code","source":"# ── 3. Temperature Scaling ────────────────────────────────────────────────\nclass TemperatureScaler(nn.Module):\n    \"\"\"Single temperature parameter T that scales all logits before sigmoid.\"\"\"\n    def __init__(self):\n        super().__init__()\n        self.temperature = nn.Parameter(torch.ones(1) * 1.5)\n\n    def forward(self, logits):\n        return logits / self.temperature\n\n\ndef fit_temperature_scaling(logits: np.ndarray, labels: np.ndarray,\n                              lr: float = 0.01, n_iter: int = 200) -> float:\n    ts_model  = TemperatureScaler()\n    optimizer = torch.optim.LBFGS([ts_model.temperature], lr=lr, max_iter=n_iter)\n    logits_t  = torch.tensor(logits, dtype=torch.float32)\n    labels_t  = torch.tensor(labels, dtype=torch.float32)\n    criterion = nn.BCEWithLogitsLoss()\n\n    def closure():\n        optimizer.zero_grad()\n        loss = criterion(ts_model(logits_t), labels_t)\n        loss.backward()\n        return loss\n\n    optimizer.step(closure)\n    T = ts_model.temperature.item()\n    return T\n\n\nT_opt = fit_temperature_scaling(raw_logits, labels)\nts_logits = raw_logits / T_opt\nts_probs  = torch.sigmoid(torch.tensor(ts_logits)).numpy()\n\nts_ece, ts_bc, ts_ba, ts_bs = expected_calibration_error(ts_probs, labels)\nprint(f'Optimal temperature T : {T_opt:.4f}')\nprint(f'Temperature-scaled ECE: {ts_ece:.5f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:47:54.070152Z","iopub.execute_input":"2026-02-20T18:47:54.070399Z","iopub.status.idle":"2026-02-20T18:48:14.644632Z","shell.execute_reply.started":"2026-02-20T18:47:54.07038Z","shell.execute_reply":"2026-02-20T18:48:14.644007Z"}},"outputs":[],"execution_count":null},{"id":"d9d63e40","cell_type":"code","source":"# ── 4. Isotonic Regression ────────────────────────────────────────────────\n# IMPORTANT: Isotonic regression is non-parametric and can perfectly memorise\n# the calibration mapping when fitted and evaluated on the same data.\n# To get an honest ECE estimate, we use 5-fold cross-validation:\n#   - Fit on 4 folds, predict on the held-out fold\n#   - Concatenate held-out predictions → compute ECE on those\n# The full-data isotonic regressor is still saved for deployment.\n\nfrom sklearn.model_selection import KFold\n\n# ── Cross-validated isotonic ECE (honest estimate) ────────────────────────\nkf = KFold(n_splits=5, shuffle=True, random_state=SEED)\nir_probs_cv = np.zeros_like(raw_probs)\n\nfor fold_idx, (train_idx, val_idx) in enumerate(kf.split(raw_probs)):\n    ir_fold = IsotonicRegression(out_of_bounds='clip')\n    ir_fold.fit(raw_probs[train_idx], labels[train_idx])\n    ir_probs_cv[val_idx] = ir_fold.predict(raw_probs[val_idx])\n\nir_ece_cv, _, _, _ = expected_calibration_error(ir_probs_cv, labels)\n\n# ── Full-data isotonic (for deployment, NOT for ECE reporting) ────────────\nir = IsotonicRegression(out_of_bounds='clip')\nir.fit(raw_probs, labels)\nir_probs = ir.predict(raw_probs)\nir_ece_naive, ir_bc, ir_ba, ir_bs = expected_calibration_error(ir_probs, labels)\n\nprint(f'Isotonic ECE (naive, same-set) : {ir_ece_naive:.5f}  ← overfitted, do NOT report')\nprint(f'Isotonic ECE (5-fold CV)       : {ir_ece_cv:.5f}  ← honest estimate')\nprint()\nprint('Note: naive ECE ≈ 0.0 is expected — isotonic memorises the mapping.')\nprint('      CV ECE is the number to compare against temperature scaling.')\n\n# Use CV ECE for fair comparison\nir_ece = ir_ece_cv\n# Use CV probs for all downstream comparison plots/tables\nir_probs_for_eval = ir_probs_cv","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:14.64564Z","iopub.execute_input":"2026-02-20T18:48:14.645935Z","iopub.status.idle":"2026-02-20T18:48:14.698352Z","shell.execute_reply.started":"2026-02-20T18:48:14.645903Z","shell.execute_reply":"2026-02-20T18:48:14.697808Z"}},"outputs":[],"execution_count":null},{"id":"ffcdc5ca","cell_type":"code","source":"# ── 4b. Overconfidence check — do extreme probabilities shrink? ────────────\nprint('=== OVERCONFIDENCE CHECK ===\\n')\n\nfor label, probs in [('Raw', raw_probs), ('Temp-scaled', ts_probs), ('Isotonic CV', ir_probs_for_eval)]:\n    n_above_99 = (probs > 0.99).sum()\n    n_below_01 = (probs < 0.01).sum()\n    n_exact_1  = (probs >= 1.0 - 1e-7).sum()\n    n_exact_0  = (probs <= 1e-7).sum()\n    print(f'{label:12s}:  p>.99 = {n_above_99:6d}  |  p<.01 = {n_below_01:6d}  |  '\n          f'p≈1.0 = {n_exact_1:5d}  |  p≈0.0 = {n_exact_0:5d}  |  '\n          f'range [{probs.min():.6f}, {probs.max():.6f}]')\n\nprint()\nif (raw_probs > 0.99).sum() > (ts_probs > 0.99).sum():\n    print('✓ Temperature scaling reduces extreme confidence — expected behaviour.')\nelse:\n    print('⚠ Temperature scaling did NOT reduce extreme confidence — check T value.')\n\n# Histogram: raw vs calibrated probability distribution\nfig, axes = plt.subplots(1, 2, figsize=(13, 4))\n\naxes[0].hist(raw_probs, bins=60, alpha=0.7, color='tab:red', label='Raw')\naxes[0].hist(ts_probs,  bins=60, alpha=0.5, color='tab:blue', label='Temp-scaled')\naxes[0].set(title='Probability distribution: Raw vs Temperature-scaled',\n            xlabel='Predicted probability', ylabel='Count')\naxes[0].legend()\n\n# Zoom into extremes (> 0.8 and < 0.2)\nextreme_mask_raw = (raw_probs > 0.8) | (raw_probs < 0.2)\nextreme_mask_ts  = (ts_probs > 0.8) | (ts_probs < 0.2)\naxes[1].hist(raw_probs[extreme_mask_raw], bins=40, alpha=0.7, color='tab:red', label='Raw (extremes)')\naxes[1].hist(ts_probs[extreme_mask_ts],   bins=40, alpha=0.5, color='tab:blue', label='Temp-scaled (extremes)')\naxes[1].set(title='Extreme tails (p<0.2 or p>0.8)',\n            xlabel='Predicted probability', ylabel='Count')\naxes[1].legend()\n\nplt.tight_layout()\nplt.savefig('/kaggle/working/overconfidence_check.png', bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:14.699296Z","iopub.execute_input":"2026-02-20T18:48:14.699917Z","iopub.status.idle":"2026-02-20T18:48:15.813938Z","shell.execute_reply.started":"2026-02-20T18:48:14.699893Z","shell.execute_reply":"2026-02-20T18:48:15.813298Z"}},"outputs":[],"execution_count":null},{"id":"da1b5f93","cell_type":"code","source":"# ── 5. Summary table & reliability diagrams ───────────────────────────────\nfrom sklearn.metrics import roc_curve\n\nsummary = pd.DataFrame([\n    {'Method': 'Raw (uncalibrated)',       'ECE': raw_ece,\n     'Brier': brier_score_loss(labels, raw_probs),\n     'AUC':   roc_auc_score(labels, raw_probs)},\n    {'Method': 'Temperature Scaling',      'ECE': ts_ece,\n     'Brier': brier_score_loss(labels, ts_probs),\n     'AUC':   roc_auc_score(labels, ts_probs)},\n    {'Method': 'Isotonic (5-fold CV)',     'ECE': ir_ece,\n     'Brier': brier_score_loss(labels, ir_probs_for_eval),\n     'AUC':   roc_auc_score(labels, ir_probs_for_eval)},\n    {'Method': 'Isotonic (naive, overfit)','ECE': ir_ece_naive,\n     'Brier': brier_score_loss(labels, ir_probs),\n     'AUC':   roc_auc_score(labels, ir_probs)},\n])\nfor col in ['ECE', 'Brier', 'AUC']:\n    summary[col] = summary[col].round(5)\n\nprint('Calibration comparison:')\nprint(summary.to_string(index=False))\n\n# ── ECE at high-confidence band (prob > 0.8) ─────────────────────────────\ndef ece_high_confidence(probs, labels, thresh=0.8):\n    \"\"\"ECE computed only for samples where predicted prob > thresh.\"\"\"\n    mask = probs > thresh\n    if mask.sum() == 0:\n        return float('nan'), 0\n    ece_val, _, _, _ = expected_calibration_error(probs[mask], labels[mask], n_bins=10)\n    return ece_val, int(mask.sum())\n\nhc_raw_ece,  hc_raw_n  = ece_high_confidence(raw_probs, labels)\nhc_ts_ece,   hc_ts_n   = ece_high_confidence(ts_probs, labels)\nhc_ir_ece,   hc_ir_n   = ece_high_confidence(ir_probs_for_eval, labels)\n\nprint(f'\\nECE at HIGH-CONFIDENCE band (p > 0.8) — critical for triage:')\nprint(f'  Raw          : ECE={hc_raw_ece:.5f}  (n={hc_raw_n:,})')\nprint(f'  Temp-scaled  : ECE={hc_ts_ece:.5f}  (n={hc_ts_n:,})')\nprint(f'  Isotonic     : ECE={hc_ir_ece:.5f}  (n={hc_ir_n:,})')\n\n# ── Reliability diagrams (3-panel) ───────────────────────────────────────\nfig, axes = plt.subplots(1, 3, figsize=(15, 5))\n\ndef plot_reliability(ax, bin_confs, bin_accs, bin_sizes, ece, title):\n    mask = bin_sizes > 0\n    ax.bar(bin_confs[mask], bin_accs[mask],\n           width=1.0/ECE_BINS, alpha=0.6, color='tab:blue',\n           align='center', label='Accuracy per bin')\n    ax.plot([0, 1], [0, 1], 'k--', linewidth=1, label='Perfect calibration')\n    gap_x = np.linspace(0, 1, 100)\n    ax.fill_between(gap_x, gap_x, alpha=0.1, color='orange', label='Calibration gap')\n    ax.set(title=f'{title}\\nECE={ece:.4f}',\n           xlabel='Mean predicted confidence', ylabel='Fraction of positives')\n    ax.set_xlim(0, 1); ax.set_ylim(0, 1)\n    ax.legend(fontsize=8)\n\nfor ax, probs, title in [(axes[0], raw_probs,         'Raw'),\n                         (axes[1], ts_probs,          'Temperature Scaled'),\n                         (axes[2], ir_probs_for_eval, 'Isotonic (CV)')]:\n    ece_i, bc_i, ba_i, bs_i = expected_calibration_error(probs, labels)\n    plot_reliability(ax, bc_i, ba_i, bs_i, ece_i, title)\n\nplt.suptitle('Reliability Diagrams', fontsize=13)\nplt.tight_layout()\nplt.savefig('/kaggle/working/reliability_diagrams.png', bbox_inches='tight')\nplt.show()\n\n# ── Overlay calibration curve: Before vs After ────────────────────────────\nfig, ax = plt.subplots(figsize=(7, 6))\n\ndef _cal_curve(probs, labels, n_bins=ECE_BINS):\n    _, confs, accs, sizes = expected_calibration_error(probs, labels, n_bins)\n    mask = sizes > 0\n    return confs[mask], accs[mask]\n\nrc, ra = _cal_curve(raw_probs, labels)\ntc, ta = _cal_curve(ts_probs, labels)\n\nax.plot([0, 1], [0, 1], 'k--', linewidth=1, label='Perfect')\nax.plot(rc, ra, 'o-', color='tab:red',  linewidth=2, markersize=5, label=f'Raw (ECE={raw_ece:.4f})')\nax.plot(tc, ta, 's-', color='tab:blue', linewidth=2, markersize=5, label=f'Temp-scaled (ECE={ts_ece:.4f})')\nax.set(title='Calibration Curve: Before vs After Temperature Scaling',\n       xlabel='Mean predicted confidence', ylabel='Fraction of positives')\nax.set_xlim(0, 1); ax.set_ylim(0, 1)\nax.legend(fontsize=10)\nax.grid(alpha=0.3)\nplt.tight_layout()\nplt.savefig('/kaggle/working/calibration_overlay.png', bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:15.814871Z","iopub.execute_input":"2026-02-20T18:48:15.815121Z","iopub.status.idle":"2026-02-20T18:48:17.072474Z","shell.execute_reply.started":"2026-02-20T18:48:15.815098Z","shell.execute_reply":"2026-02-20T18:48:17.071874Z"}},"outputs":[],"execution_count":null},{"id":"4b6540a9","cell_type":"code","source":"# ── 6. Choose best calibration method & re-optimise threshold ─────────────\n# Use Temperature Scaling as the deployed method:\n#   - Simple (1 parameter), generalises well, preserves AUC\n#   - Isotonic can overfit on smaller val sets\n# Change to 'isotonic' if ir_ece is substantially better\n\nCALIBRATION_METHOD = 'temperature'   # options: 'none', 'temperature', 'isotonic'\n\nif CALIBRATION_METHOD == 'temperature':\n    cal_probs = ts_probs\n    best_ece  = ts_ece\nelif CALIBRATION_METHOD == 'isotonic':\n    cal_probs = ir_probs\n    best_ece  = ir_ece\nelse:\n    cal_probs = raw_probs\n    best_ece  = raw_ece\n\n# ── Threshold re-optimisation after calibration ──────────────────────────\n# Calibration shifts the probability distribution, so the old threshold\n# (Youden-optimised on raw probs) is no longer optimal.  We recompute.\nfpr, tpr, thresholds = roc_curve(labels, cal_probs)\nj_scores = tpr - fpr\nbest_idx = np.argmax(j_scores)\nCAL_THRESHOLD = float(thresholds[best_idx])\n\nprint(f'Threshold comparison:')\nprint(f'  Original (raw, from NB03) : {BASE_THRESHOLD:.4f}')\nprint(f'  Re-optimised (calibrated) : {CAL_THRESHOLD:.4f}')\nprint(f'  Shift                     : {CAL_THRESHOLD - BASE_THRESHOLD:+.4f}')\nprint()\n\n# Show the impact on sensitivity/specificity\ncal_preds_old = (cal_probs >= BASE_THRESHOLD).astype(int)\ncal_preds_new = (cal_probs >= CAL_THRESHOLD).astype(int)\ntp_old = ((cal_preds_old == 1) & (labels == 1)).sum()\ntp_new = ((cal_preds_new == 1) & (labels == 1)).sum()\nfp_old = ((cal_preds_old == 1) & (labels == 0)).sum()\nfp_new = ((cal_preds_new == 1) & (labels == 0)).sum()\nfn_old = ((cal_preds_old == 0) & (labels == 1)).sum()\nfn_new = ((cal_preds_new == 0) & (labels == 1)).sum()\ntn_old = ((cal_preds_old == 0) & (labels == 0)).sum()\ntn_new = ((cal_preds_new == 0) & (labels == 0)).sum()\nsens_old = tp_old / max(tp_old + fn_old, 1)\nsens_new = tp_new / max(tp_new + fn_new, 1)\nspec_old = tn_old / max(tn_old + fp_old, 1)\nspec_new = tn_new / max(tn_new + fp_new, 1)\n\nprint(f'  Using old threshold → Sens={sens_old:.4f}  Spec={spec_old:.4f}')\nprint(f'  Using new threshold → Sens={sens_new:.4f}  Spec={spec_new:.4f}')\n\n# ── Sensitivity-constrained threshold (clinical objective) ───────────────\n# For hemorrhage detection, false negatives are more costly than false positives.\n# Clinical approach: fix sensitivity ≥ target, then maximise specificity.\nSENS_TARGETS = [0.95, 0.97, 0.99]\nprint(f'\\nSensitivity-constrained thresholds (on calibrated probs):')\nprint(f'  {\"Target Sens\":>12s}  {\"Threshold\":>10s}  {\"Actual Sens\":>11s}  {\"Specificity\":>11s}')\nfor sens_target in SENS_TARGETS:\n    # Find highest threshold where tpr >= target\n    valid = tpr >= sens_target\n    if valid.any():\n        # Among valid, pick the one with lowest fpr (highest specificity)\n        idx = np.where(valid)[0]\n        best = idx[np.argmin(fpr[idx])]\n        thr = float(thresholds[best])\n        preds = (cal_probs >= thr).astype(int)\n        tp_s = ((preds == 1) & (labels == 1)).sum()\n        fn_s = ((preds == 0) & (labels == 1)).sum()\n        tn_s = ((preds == 0) & (labels == 0)).sum()\n        fp_s = ((preds == 1) & (labels == 0)).sum()\n        s = tp_s / max(tp_s + fn_s, 1)\n        sp = tn_s / max(tn_s + fp_s, 1)\n        print(f'  {sens_target:>11.0%}  {thr:>10.4f}  {s:>10.4f}  {sp:>10.4f}')\n    else:\n        print(f'  {sens_target:>11.0%}  {\"N/A\":>10s}  {\"—\":>11s}  {\"—\":>11s}')\n\n# ── 3-band triage system ──────────────────────────────────────────────────\n#  HIGH   : cal_prob >= 0.75  → haemorrhage detected, URGENT review\n#  MEDIUM : 0.35 <= cal_prob < 0.75 → uncertain, standard review\n#  LOW    : cal_prob < 0.35   → likely normal, routine workflow\n\nHIGH_THRESHOLD   = 0.75\nLOW_THRESHOLD    = 0.35\n\ndef assign_band(p):\n    if p >= HIGH_THRESHOLD:\n        return 'HIGH'\n    elif p >= LOW_THRESHOLD:\n        return 'MEDIUM'\n    else:\n        return 'LOW'\n\nbands = np.vectorize(assign_band)(cal_probs)\n\nband_df = pd.DataFrame({\n    'label'   : labels,\n    'cal_prob': cal_probs,\n    'band'    : bands,\n})\n\nprint(f'\\nBand distribution:')\nprint(band_df['band'].value_counts().to_frame())\nprint(f'\\n{CALIBRATION_METHOD.upper()} ECE: {best_ece:.5f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:17.073443Z","iopub.execute_input":"2026-02-20T18:48:17.073768Z","iopub.status.idle":"2026-02-20T18:48:17.102682Z","shell.execute_reply.started":"2026-02-20T18:48:17.073738Z","shell.execute_reply":"2026-02-20T18:48:17.102118Z"}},"outputs":[],"execution_count":null},{"id":"158bfea4","cell_type":"code","source":"# ── 7. Per-band error rates ───────────────────────────────────────────────\nband_stats = []\nfor band_name in ['HIGH', 'MEDIUM', 'LOW']:\n    sub = band_df[band_df['band'] == band_name]\n    if len(sub) == 0: continue\n    n_pos      = int(sub['label'].sum())\n    n_neg      = int((sub['label'] == 0).sum())\n    n_total    = len(sub)\n    pos_rate   = n_pos / n_total\n\n    # For HIGH band: FN rate (fraction of actual positives predicted in a lower band)\n    # For LOW  band: FP rate (fraction of actual negatives predicted in a higher band)\n    # Here we just report positive rate and fraction of cases in each band.\n    band_stats.append({\n        'Band'         : band_name,\n        'N cases'      : n_total,\n        '% of val set' : round(n_total / len(band_df) * 100, 1),\n        'Positive rate': round(pos_rate, 4),\n        'N positive'   : n_pos,\n        'N negative'   : n_neg,\n    })\n\nband_summary = pd.DataFrame(band_stats)\nprint('\\nPer-band statistics:')\nprint(band_summary.to_string(index=False))\nprint()\nprint('Interpretation:')\nprint(f'  HIGH band positive rate ≈ PPV for urgent review triage')\nprint(f'  LOW  band positive rate ≈ miss rate for de-prioritised cases')\nprint(f'  MEDIUM band = requires manual review — maximise routing here for safety')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:17.103558Z","iopub.execute_input":"2026-02-20T18:48:17.103804Z","iopub.status.idle":"2026-02-20T18:48:17.121089Z","shell.execute_reply.started":"2026-02-20T18:48:17.103782Z","shell.execute_reply":"2026-02-20T18:48:17.120451Z"}},"outputs":[],"execution_count":null},{"id":"02ae1411","cell_type":"code","source":"# ── 8. Visualise confidence band distribution ─────────────────────────────\nfig, axes = plt.subplots(1, 2, figsize=(12, 5))\n\nband_colors = {'HIGH': '#e74c3c', 'MEDIUM': '#f39c12', 'LOW': '#2ecc71'}\n\n# Histogram of calibrated probs, coloured by band\nfor band_name, color in band_colors.items():\n    sub = band_df[band_df['band'] == band_name]['cal_prob']\n    axes[0].hist(sub, bins=40, alpha=0.7, color=color,\n                 label=f'{band_name} (n={len(sub)})')\naxes[0].axvline(HIGH_THRESHOLD, color='red',   linestyle='--', linewidth=1)\naxes[0].axvline(LOW_THRESHOLD,  color='green', linestyle='--', linewidth=1)\naxes[0].set(title='Calibrated probability distribution by band',\n            xlabel='Calibrated probability', ylabel='Count')\naxes[0].legend()\n\n# Positive rate per band\naxes[1].bar(band_summary['Band'], band_summary['Positive rate'],\n            color=[band_colors[b] for b in band_summary['Band']])\nfor i, row in band_summary.iterrows():\n    axes[1].text(i, row['Positive rate'] + 0.01,\n                 f'{row[\"Positive rate\"]*100:.1f}%',\n                 ha='center', fontsize=10, fontweight='bold')\naxes[1].set(title='Hemorrhage positive rate per confidence band',\n            ylabel='Positive rate', ylim=(0, 1.1))\n\nplt.suptitle(f'Confidence Band Analysis ({CALIBRATION_METHOD.title()} calibration)', fontsize=12)\nplt.tight_layout()\nplt.savefig('/kaggle/working/confidence_bands.png', bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:17.123326Z","iopub.execute_input":"2026-02-20T18:48:17.123606Z","iopub.status.idle":"2026-02-20T18:48:17.876155Z","shell.execute_reply.started":"2026-02-20T18:48:17.123585Z","shell.execute_reply":"2026-02-20T18:48:17.875556Z"}},"outputs":[],"execution_count":null},{"id":"bdfdb5a6","cell_type":"code","source":"# ── 9. Save calibration parameters ───────────────────────────────────────\nimport joblib\n\nraw_brier = brier_score_loss(labels, raw_probs)\ncal_brier = brier_score_loss(labels, cal_probs)\n\ncalibration_params = {\n    'method'              : CALIBRATION_METHOD,\n    'temperature'         : float(T_opt),\n    'base_threshold'      : float(BASE_THRESHOLD),\n    'calibrated_threshold': float(CAL_THRESHOLD),\n    'high_threshold'      : float(HIGH_THRESHOLD),\n    'low_threshold'       : float(LOW_THRESHOLD),\n    'raw_ece'             : float(raw_ece),\n    'cal_ece'             : float(best_ece),\n    'raw_brier'           : float(raw_brier),\n    'cal_brier'           : float(cal_brier),\n    'hc_ece_raw'          : float(hc_raw_ece),\n    'hc_ece_cal'          : float(hc_ts_ece if CALIBRATION_METHOD == 'temperature' else hc_ir_ece),\n}\n\nwith open('/kaggle/working/calibration_params.json', 'w') as f:\n    json.dump(calibration_params, f, indent=2)\n\n# Save isotonic regressor (useful as fallback)\njoblib.dump(ir, '/kaggle/working/isotonic_regressor.pkl')\n\nprint('Saved:')\nprint('  /kaggle/working/calibration_params.json')\nprint('  /kaggle/working/isotonic_regressor.pkl')\nprint('  /kaggle/working/reliability_diagrams.png')\nprint('  /kaggle/working/calibration_overlay.png')\nprint('  /kaggle/working/overconfidence_check.png')\nprint('  /kaggle/working/confidence_bands.png')\nprint()\nprint(json.dumps(calibration_params, indent=2))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:17.877064Z","iopub.execute_input":"2026-02-20T18:48:17.877339Z","iopub.status.idle":"2026-02-20T18:48:17.887228Z","shell.execute_reply.started":"2026-02-20T18:48:17.877305Z","shell.execute_reply":"2026-02-20T18:48:17.88669Z"}},"outputs":[],"execution_count":null},{"id":"967d0511","cell_type":"code","source":"# ── HEALTH CHECK — automated output validation ────────────────────────────\n\nerrors = []\n\n# Check calibration params\ncal_hc_path = '/kaggle/working/calibration_params.json'\nif not os.path.exists(cal_hc_path):\n    errors.append('calibration_params.json is MISSING')\nelse:\n    with open(cal_hc_path) as f:\n        cal_hc = json.load(f)\n    if cal_hc.get('temperature', 0) <= 0:\n        errors.append(f'Invalid temperature: {cal_hc.get(\"temperature\")}')\n    if cal_hc.get('cal_ece', 1) > 0.20:\n        errors.append(f'Calibrated ECE={cal_hc[\"cal_ece\"]:.4f} is high — calibration may have failed')\n    if 'calibrated_threshold' not in cal_hc:\n        errors.append('calibrated_threshold missing from params — threshold was not re-optimised')\n\n# Check isotonic regressor\nif not os.path.exists('/kaggle/working/isotonic_regressor.pkl'):\n    errors.append('isotonic_regressor.pkl is MISSING')\n\n# Check plots\nfor plot in ['reliability_diagrams.png', 'confidence_bands.png',\n             'calibration_overlay.png', 'overconfidence_check.png']:\n    if not os.path.exists(f'/kaggle/working/{plot}'):\n        errors.append(f'Missing plot: {plot}')\n\n# Caveat reminder\nprint('─' * 60)\nprint('CALIBRATION CAVEATS:')\nprint('  1. Temperature and isotonic were fitted AND evaluated on the')\nprint('     SAME validation set → post-calibration ECE is slightly')\nprint('     optimistic.')\nprint('  2. Calibration is at SLICE level. Patient-level aggregation')\nprint('     (e.g. max-pooling) changes the probability distribution')\nprint('     and may require re-calibration.')\nprint('─' * 60)\n\nhealth = {\n    'notebook'            : '06_calibration',\n    'status'              : 'PASS' if not errors else 'FAIL',\n    'errors'              : errors,\n    'raw_ece'             : round(float(raw_ece), 5),\n    'cal_ece'             : round(float(best_ece), 5),\n    'raw_brier'           : round(float(raw_brier), 5),\n    'cal_brier'           : round(float(cal_brier), 5),\n    'temperature'         : round(float(T_opt), 4),\n    'base_threshold'      : round(float(BASE_THRESHOLD), 4),\n    'calibrated_threshold': round(float(CAL_THRESHOLD), 4),\n    'method'              : CALIBRATION_METHOD,\n    'hc_ece_raw'          : round(float(hc_raw_ece), 5),\n    'hc_ece_cal'          : round(float(hc_ts_ece if CALIBRATION_METHOD == 'temperature' else hc_ir_ece), 5),\n    'caveats'             : ['Fitted and evaluated on same val set',\n                             'Slice-level calibration only'],\n}\n\nwith open('/kaggle/working/health_check_nb06.json', 'w') as f:\n    json.dump(health, f, indent=2)\n\nif errors:\n    print('❌ HEALTH CHECK FAILED:')\n    for e in errors:\n        print(f'   • {e}')\nelse:\n    print('✅ HEALTH CHECK PASSED')\n    print(f'   Method              : {CALIBRATION_METHOD}')\n    print(f'   Temperature         : {T_opt:.4f}')\n    print(f'   Raw ECE / Brier     : {raw_ece:.5f} / {raw_brier:.5f}')\n    print(f'   Cal ECE / Brier     : {best_ece:.5f} / {cal_brier:.5f}')\n    print(f'   HC-band ECE (p>.8)  : raw={hc_raw_ece:.5f}  cal={hc_ts_ece:.5f}')\n    print(f'   Original threshold  : {BASE_THRESHOLD:.4f}')\n    print(f'   Calibrated threshold: {CAL_THRESHOLD:.4f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-20T18:48:17.888728Z","iopub.execute_input":"2026-02-20T18:48:17.888984Z","iopub.status.idle":"2026-02-20T18:48:17.906514Z","shell.execute_reply.started":"2026-02-20T18:48:17.888953Z","shell.execute_reply":"2026-02-20T18:48:17.905853Z"}},"outputs":[],"execution_count":null}]}