{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "39244d24",
   "metadata": {},
   "source": [
    "# Classification Metrics for Credit\n",
    "\n",
    "Companion notebook to [Data Science Module 4](https://rodoslay.com/data-science/04-classification-metrics) at **rodoslay.com**.\n",
    "\n",
    "Every metric built from scratch before touching `sklearn`: ROC and PR curves by hand, Gini three\n",
    "ways (they agree — proven, not asserted), KS with its optimal cutoff, reliability diagrams with\n",
    "isotonic recalibration, and the **base-rate experiment**: one model, three prevalence levels,\n",
    "a table of which metrics moved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd65c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.isotonic import IsotonicRegression\n",
    "\n",
    "rng = np.random.default_rng(4)\n",
    "plt.rcParams[\"figure.figsize\"] = (10, 4.5)\n",
    "plt.rcParams[\"axes.grid\"] = True\n",
    "plt.rcParams[\"grid.alpha\"] = 0.3\n",
    "\n",
    "def make_scores(n, base_rate, sep=1.5):\n",
    "    \"\"\"Goods ~ N(0,1), bads ~ N(sep,1); returns scores and labels.\"\"\"\n",
    "    n_bad = int(n * base_rate)\n",
    "    y = np.r_[np.zeros(n - n_bad, int), np.ones(n_bad, int)]\n",
    "    s = np.r_[rng.normal(0, 1, n - n_bad), rng.normal(sep, 1, n_bad)]\n",
    "    idx = rng.permutation(n)\n",
    "    return s[idx], y[idx]\n",
    "\n",
    "scores, y = make_scores(50_000, base_rate=0.05)\n",
    "print(f\"n = {len(y):,}, base rate = {y.mean():.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc20d0a4",
   "metadata": {},
   "source": [
    "## 1. ROC, by hand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86502c40",
   "metadata": {},
   "outputs": [],
   "source": [
    "order = np.argsort(-scores)                     # descending score = most suspicious first\n",
    "y_sorted = y[order]\n",
    "tpr = np.cumsum(y_sorted) / y.sum()             # bads caught so far\n",
    "fpr = np.cumsum(1 - y_sorted) / (1 - y).sum()   # goods falsely flagged so far\n",
    "auc = np.trapezoid(tpr, fpr)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5.5, 5))\n",
    "ax.plot(fpr, tpr, color=\"#1f4d3a\"); ax.plot([0, 1], [0, 1], \"k--\", lw=0.8)\n",
    "ax.set_xlabel(\"FPR\"); ax.set_ylabel(\"TPR\"); ax.set_title(f\"ROC by hand — AUC = {auc:.4f}\")\n",
    "plt.show()\n",
    "\n",
    "# check the probabilistic meaning directly: P(random bad scores above random good)\n",
    "bads, goods = scores[y == 1], scores[y == 0]\n",
    "pairs = rng.integers(0, [len(bads), len(goods)], size=(200_000, 2))\n",
    "print(f\"P(bad > good) by sampling: {(bads[pairs[:,0]] > goods[pairs[:,1]]).mean():.4f}  (= AUC)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a9e5820",
   "metadata": {},
   "source": [
    "## 2. Gini three ways"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e252d6b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# way 1: from AUC\n",
    "gini_1 = 2 * auc - 1\n",
    "\n",
    "# way 2: CAP / accuracy ratio — cumulative bads captured vs share of population flagged\n",
    "pop = np.arange(1, len(y) + 1) / len(y)\n",
    "cap = np.cumsum(y_sorted) / y.sum()\n",
    "a_model = np.trapezoid(cap, pop) - 0.5\n",
    "a_perfect = (1 - y.mean() / 2) - 0.5\n",
    "gini_2 = a_model / a_perfect\n",
    "\n",
    "# way 3: Somers' D via concordance sampling\n",
    "b = bads[rng.integers(0, len(bads), 300_000)]\n",
    "g = goods[rng.integers(0, len(goods), 300_000)]\n",
    "gini_3 = (b > g).mean() - (b < g).mean()\n",
    "\n",
    "print(f\"Gini via AUC:        {gini_1:.4f}\")\n",
    "print(f\"Gini via CAP curve:  {gini_2:.4f}\")\n",
    "print(f\"Gini via Somers' D:  {gini_3:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32f2c6f6",
   "metadata": {},
   "source": [
    "Three constructions, one number (up to sampling noise) — quote whichever your audience expects.\n",
    "\n",
    "## 3. KS and its cutoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a07912a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = np.linspace(scores.min(), scores.max(), 400)\n",
    "cdf_good = np.searchsorted(np.sort(goods), grid) / len(goods)\n",
    "cdf_bad = np.searchsorted(np.sort(bads), grid) / len(bads)\n",
    "ks = (cdf_good - cdf_bad).max()\n",
    "ks_at = grid[np.argmax(cdf_good - cdf_bad)]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(grid, cdf_good, label=\"goods CDF\", color=\"#1f4d3a\")\n",
    "ax.plot(grid, cdf_bad, label=\"bads CDF\", color=\"#b3541e\")\n",
    "ax.vlines(ks_at, cdf_bad[np.argmax(cdf_good - cdf_bad)], cdf_good[np.argmax(cdf_good - cdf_bad)],\n",
    "          color=\"k\", ls=\"--\", label=f\"KS = {ks:.3f}\")\n",
    "ax.legend(); ax.set_title(\"KS: the best single cutoff's separation\"); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e5b8dc9",
   "metadata": {},
   "source": [
    "## 4. The base-rate experiment\n",
    "\n",
    "Same separation (the same *model*), three prevalence levels. Threshold fixed at the KS-optimal\n",
    "point each time. Watch the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e09e80",
   "metadata": {},
   "outputs": [],
   "source": [
    "rows = []\n",
    "for br in [0.50, 0.10, 0.02]:\n",
    "    s, yy = make_scores(60_000, base_rate=br)\n",
    "    gd, bd = s[yy == 0], s[yy == 1]\n",
    "    # ranking metrics\n",
    "    o = np.argsort(-s); ys = yy[o]\n",
    "    tp_r = np.cumsum(ys) / yy.sum(); fp_r = np.cumsum(1 - ys) / (1 - yy).sum()\n",
    "    auc_ = np.trapezoid(tp_r, fp_r)\n",
    "    gr = np.linspace(s.min(), s.max(), 400)\n",
    "    ks_ = (np.searchsorted(np.sort(gd), gr) / len(gd) - np.searchsorted(np.sort(bd), gr) / len(bd)).max()\n",
    "    thr = gr[np.argmax(np.searchsorted(np.sort(gd), gr) / len(gd) - np.searchsorted(np.sort(bd), gr) / len(bd))]\n",
    "    # threshold metrics\n",
    "    flag = s > thr\n",
    "    tp = (flag & (yy == 1)).sum(); fp = (flag & (yy == 0)).sum()\n",
    "    fn = (~flag & (yy == 1)).sum(); tn = (~flag & (yy == 0)).sum()\n",
    "    acc = (tp + tn) / len(yy); prec = tp / max(tp + fp, 1); rec = tp / max(tp + fn, 1)\n",
    "    naive = max(br, 1 - br)\n",
    "    rows.append([f\"{br:.0%}\", f\"{auc_:.3f}\", f\"{ks_:.3f}\", f\"{acc:.1%}\", f\"{naive:.1%}\", f\"{prec:.1%}\", f\"{rec:.1%}\"])\n",
    "\n",
    "pd.DataFrame(rows, columns=[\"base rate\", \"AUC\", \"KS\", \"accuracy\", \"accuracy of 'approve all'\", \"precision\", \"recall\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e85e0b9d",
   "metadata": {},
   "source": [
    "The table is the module's argument in numbers: **AUC and KS barely move** (they're base-rate\n",
    "free); **accuracy rises as the model becomes less useful**, converging to the approve-everyone\n",
    "baseline; **precision collapses** — at 2% prevalence, most flags are false alarms even though the\n",
    "model is unchanged. Decide with the metric that feels the base rate.\n",
    "\n",
    "## 5. Calibration and isotonic recalibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b47d9986",
   "metadata": {},
   "outputs": [],
   "source": [
    "# miscalibrate on purpose: train-style scores -> overconfident \"probabilities\"\n",
    "s, yy = make_scores(40_000, base_rate=0.05)\n",
    "p_bad_prior = 1 / (1 + np.exp(-(1.9 * s - 2.0)))      # a badly scaled sigmoid\n",
    "\n",
    "def reliability(p, y, bins=10):\n",
    "    q = np.quantile(p, np.linspace(0, 1, bins + 1)); q[0], q[-1] = 0, 1\n",
    "    idx = np.digitize(p, q[1:-1])\n",
    "    return np.array([[p[idx == b].mean(), y[idx == b].mean()] for b in range(bins) if (idx == b).any()])\n",
    "\n",
    "half = len(yy) // 2                                    # fit map on first half, evaluate on second\n",
    "iso = IsotonicRegression(out_of_bounds=\"clip\").fit(p_bad_prior[:half], yy[:half])\n",
    "p_cal = iso.predict(p_bad_prior[half:])\n",
    "\n",
    "brier_raw = np.mean((p_bad_prior[half:] - yy[half:]) ** 2)\n",
    "brier_cal = np.mean((p_cal - yy[half:]) ** 2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "for probs, label, color in [(p_bad_prior[half:], f\"raw (Brier {brier_raw:.4f})\", \"#b3541e\"),\n",
    "                            (p_cal, f\"isotonic (Brier {brier_cal:.4f})\", \"#1f4d3a\")]:\n",
    "    r = reliability(probs, yy[half:])\n",
    "    ax.plot(r[:, 0], r[:, 1], \"o-\", label=label, color=color)\n",
    "lim = max(p_bad_prior[half:].max(), 0.5)\n",
    "ax.plot([0, lim], [0, lim], \"k--\", lw=0.8)\n",
    "ax.set_xlabel(\"predicted PD\"); ax.set_ylabel(\"realized\"); ax.legend()\n",
    "ax.set_title(\"Recalibration fixes the level; the ranking (and Gini) is untouched\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5700d883",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. **PR vs ROC.** Plot both curves for the 2% base-rate dataset. Which one would have warned the\n",
    "   collections team about the false-alarm rate?\n",
    "2. **The pricing shoot-out.** Model A: separation 1.6 but overconfident probabilities. Model B:\n",
    "   separation 1.4, perfectly calibrated. Price loans at rate = 5% + PD×LGD/(1−PD) with LGD = 45%,\n",
    "   fund at 3%, and compute realized profit under each. Which wins, and what does that say about\n",
    "   Gini-only model selection?\n",
    "3. **Confidence intervals.** Bootstrap the AUC (resample rows 1,000 times). How wide is the 95%\n",
    "   interval at n = 5,000 — and is your favorite \"improvement\" of +0.005 AUC even real?\n",
    "4. **Segment slices.** Add a thin-file flag and compute Gini per segment. Construct a case where\n",
    "   portfolio Gini rises while thin-file Gini falls."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}