{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "interpretableml.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# interpretable machine learning (interpretable ML) \u2014 Python demo\n\nNumerical companion to the entry [interpretable machine learning (interpretable ML)](https://dictionaryofml.org/terms/interpretableml.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nBacks the entry's two empirical claims on a synthetic trainset: restricting the hypothesis space to interpretable (sparse) hypotheses acts as regularization, and the resulting interpretable hypothesis achieves accuracy comparable to an opaque (dense) one. Also verifies that the entry's triage decision tree is an explicit, traceable function. Self-contained (numpy only), fixed seeds.\n\nRequires NumPy only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/interpretableml.py`](https://dictionaryofml.org/terms/interpretableml.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"interpretableml.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ninterpretableml.py \u2014 numerical companion to the glossary entry\n'interpretable machine learning (interpretable ML)'.\n\nPurpose\n-------\nBacks the entry's two empirical claims on a synthetic trainset:\nrestricting the hypothesis space to interpretable (sparse) hypotheses\nacts as regularization, and the resulting interpretable hypothesis\nachieves accuracy comparable to an opaque (dense) one.  Also verifies\nthat the entry's triage decision tree is an explicit, traceable\nfunction.  Self-contained (numpy only), fixed seeds.\n\nBlocks\n------\n[B-reg]   Regularization by pruning the model: on m = 30 training\n          points with d = 20 features (3 informative), the dense\n          linear-regression hypothesis attains a lower error on the\n          trainset than the sparse one (it also fits the 17 noise\n          features), but its TEST error is higher \u2014 restricting the\n          hypothesis space to sparse, interpretable hypotheses lowers\n          the test error.\n[B-acc]   Accuracy parity: the sparse hypothesis uses at most 6 of the\n          20 features (a human can read its few terms) and its test\n          error is at most that of the dense hypothesis.\n[B-tree]  The entry's triage tree is its own explanation: the tree\n          prediction (temperature > 39 -> urgent; else heart rate\n          > 120 -> urgent; else routine) coincides with the explicit\n          two-rule lookup on a grid of 10000 (temperature, heart rate)\n          pairs, and it partitions the possible (temperature, heart\n          rate) pairs into exactly three axis-aligned cells.\n[B-effect] The two effects an interpretable hypothesis makes\n          traceable, both exact for a linear model fitted by the normal\n          equations: moving feature j by Delta moves the prediction by\n          w_j * Delta, and perturbing training label r by delta shifts a\n          prediction by delta * x^T (X^T X)^{-1} x^(r) -- computable\n          without refitting, which refitting then confirms.\n[B-lime]  LIME idea: a proximity-weighted linear fit approximates a\n          nonlinear learned hypothesis near the data point x0 = 0.8\n          (max error < 0.3 within |x - x0| < 0.3) but not globally\n          (error more than 4x larger far from x0).  In d = 5 features, the local\n          linear surrogate of a hypothesis that depends only on the\n          first two features has near-zero coefficients for the other\n          three: the surrogate depends on few features.\n[B-eerm]  EERM idea (Zhang et al., 2024): the penalty charges the\n          departure of the learned hypothesis from the predictions\n          one user supplies.  Those predictions are not elicited one\n          data point at a time here; they are the predictions of a\n          simpler proxy model the user considers interpretable, a\n          linear-regression line.  Adding\n          lambda * mean((h(x) - userpred(x))^2) to the training of a\n          degree-10 polynomial pulls the learned\n          hypothesis toward that line; the EERM fit has (i) smaller\n          departure from the user predictions and (ii) smaller test\n          error than the fit without the penalty, at a slightly\n          larger error on the trainset.\n\nThe sparse hypothesis is the relaxed Lasso under a design with\nX^T X = I: the Lasso there has the closed form of soft-thresholding\nthe linear-regression coefficients (Hastie et al., 2009,\nSect. 3.4.2), it is used to select the features, and linear\nregression is refit on the selected support; no iterative solver is\nneeded.\n\nOutputs\n-------\ninterpretableml_lime.csv     : x, h (learned hypothesis), g (local\n                               linear surrogate; nan outside the\n                               window |x - x0| <= 0.6) for the LIME\n                               figure.\ninterpretableml_eerm.csv     : x, hopaque (unregularized polynomial),\n                               hsurr (the user predictions, from the\n                               proxy line), heerm\n                               (EERM-regularized polynomial) on a\n                               grid, for the EERM figure.\ninterpretableml_eermdata.csv : x, y training points of the EERM\n                               figure.\n(The decision-tree figure is schematic TikZ and needs no data.)\n\"\"\"\n\nimport numpy as np                  # the only dependency\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []                         # collects (check name, pass/fail) pairs\n\n\ndef check(name, ok):                # records and prints one verification\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\nrng = np.random.default_rng(4)\n\n# synthetic tabular regression: d = 20 features, only 3 informative\nd, m_train = 20, 30\nw_true = np.zeros(d)\nw_true[:3] = (3.0, -2.0, 1.5)       # the interpretable ground truth\nsigma = 0.5                          # label noise level\n\n# design with X^T X = I, so linear regression and Lasso are closed form\nX, _ = np.linalg.qr(rng.standard_normal((m_train, d)))\ny = X @ w_true + sigma * rng.standard_normal(m_train)\n\nw_dense = X.T @ y                    # dense linear-regression hypothesis\nlam = 1.0                            # Lasso penalty weight\n# Lasso under X^T X = I = soft-thresholding of w_dense; the Lasso\n# here SELECTS the features, and the sparse hypothesis refits linear\n# regression on the selected support (relaxed Lasso) \u2014 with X^T X = I\n# that refit just keeps the selected entries of w_dense\nsupport = np.abs(w_dense) > lam\nw_sparse = np.where(support, w_dense, 0.0)\n\n# errors on the trainset: dense fits the noise features, so it fits better\ntrain_dense = np.mean((y - X @ w_dense) ** 2)\ntrain_sparse = np.mean((y - X @ w_sparse) ** 2)\n# expected test error for x ~ N(0, I): ||w - w_true||^2 + sigma^2\ntest_dense = np.sum((w_dense - w_true) ** 2) + sigma ** 2\ntest_sparse = np.sum((w_sparse - w_true) ** 2) + sigma ** 2"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-reg]** Regularization by pruning the model: on m = 30 training points with d = 20 features (3 informative), the dense linear-regression hypothesis attains a lower error on the trainset than the sparse one (it also fits the 17 noise features), but its TEST error is higher \u2014 restricting the hypothesis space to sparse, interpretable hypotheses lowers the test error."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "check(\"[B-reg]   dense fits training data better but generalizes worse\",\n      train_dense < train_sparse and test_sparse < test_dense)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-acc]** Accuracy parity: the sparse hypothesis uses at most 6 of the 20 features (a human can read its few terms) and its test error is at most that of the dense hypothesis."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "n_used = int(np.count_nonzero(w_sparse))\ncheck(\"[B-acc]   sparse hypothesis uses few features at comparable \"\n      f\"test error ({n_used} of {d})\",\n      n_used <= 6 and test_sparse <= test_dense)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-tree]** The entry's triage tree is its own explanation: the tree prediction (temperature > 39 -> urgent; else heart rate > 120 -> urgent; else routine) coincides with the explicit two-rule lookup on a grid of 10000 (temperature, heart rate) pairs, and it partitions the possible (temperature, heart rate) pairs into exactly three axis-aligned cells."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def tree(temp, hr):                  # the entry's triage decision tree\n    if temp > 39.0:\n        return \"urgent\"\n    if hr > 120.0:\n        return \"urgent\"\n    return \"routine\"\n\n\ndef two_rules(temp, hr):             # the explicit lookup a human reads off\n    return \"urgent\" if (temp > 39.0 or hr > 120.0) else \"routine\"\n\n\ntemps = np.linspace(35.0, 42.0, 100)\nhrs = np.linspace(50.0, 180.0, 100)\nagree = all(tree(t, h) == two_rules(t, h) for t in temps for h in hrs)\n# the tree partitions the plane into 3 axis-aligned cells:\n# {temp > 39}, {temp <= 39, hr > 120}, {temp <= 39, hr <= 120}\ncells = {(t > 39.0, t <= 39.0 and h > 120.0) for t in temps for h in hrs}\ncheck(\"[B-tree]  tree prediction equals the two-rule lookup on a \"\n      \"100 x 100 grid; 3 cells\", agree and len(cells) == 3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-effect]** The two effects an interpretable hypothesis makes traceable, both exact for a linear model fitted by the normal equations: moving feature j by Delta moves the prediction by w_j * Delta, and perturbing training label r by delta shifts a prediction by delta * x^T (X^T X)^{-1} x^(r) -- computable without refitting, which refitting then confirms."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# The two effects the entry claims a user of a linear model can work out:\n# that of a change to one feature, and that of a change to one training\n# label.\n# Both are exact here, not approximate, which is what \"traceable\" means.\nrng_s = np.random.default_rng(11)\nm_s, d_s = 12, 3\nXs = rng_s.standard_normal((m_s, d_s))\nys = Xs @ np.array([1.5, -0.8, 0.3]) + 0.2 * rng_s.standard_normal(m_s)\n\nG = Xs.T @ Xs                                      # the normal equations\nw_hat = np.linalg.solve(G, Xs.T @ ys)\n\n# (i) feature change: moving feature j by Delta moves the prediction by\n# exactly w_j * Delta, whatever the other features hold\nx0 = rng_s.standard_normal(d_s)\nj, delta_x = 1, 0.7\nx1 = x0.copy()\nx1[j] += delta_x\nmoved = float(x1 @ w_hat - x0 @ w_hat)\nprint(f\"  feature {j} moved by {delta_x}: prediction moves {moved:+.4f}, \"\n      f\"weight times step {w_hat[j] * delta_x:+.4f}\")\ncheck(\"[B-effect] a feature change moves the prediction by weight times step\",\n      abs(moved - w_hat[j] * delta_x) < 1e-12)\n\n# (ii) label change: the normal equations make the map from training set to\n# hypothesis LINEAR in the labels, so perturbing label r by delta shifts a\n# prediction by exactly delta * x^T (X^T X)^{-1} x^(r).  The user can read\n# that off without refitting; refitting here only confirms it.\nr, delta_y = 4, 0.5\npredicted = float(delta_y * x0 @ np.linalg.solve(G, Xs[r]))\nys_pert = ys.copy()\nys_pert[r] += delta_y\nw_pert = np.linalg.solve(G, Xs.T @ ys_pert)\nactual = float(x0 @ w_pert - x0 @ w_hat)\nprint(f\"  label {r} perturbed by {delta_y}: prediction shifts {actual:+.4f}, \"\n      f\"normal equations predict {predicted:+.4f}\")\ncheck(\"[B-effect] a perturbed training label shifts the prediction by the \"\n      \"amount the normal equations predict\", abs(actual - predicted) < 1e-12)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-lime]** LIME idea: a proximity-weighted linear fit approximates a nonlinear learned hypothesis near the data point x0 = 0.8 (max error < 0.3 within |x - x0| < 0.3) but not globally (error more than 4x larger far from x0). In d = 5 features, the local linear surrogate of a hypothesis that depends only on the first two features has near-zero coefficients for the other three: the surrogate depends on few features."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# LIME idea: a proximity-weighted linear fit approximates the learned\n# hypothesis near the data point x0, using few features.\nrng_l = np.random.default_rng(5)\n\n\ndef h_hat(x):                        # stand-in for an opaque learned hypothesis\n    return np.sin(3.0 * x) + 0.5 * x ** 2\n\n\nx0 = 0.8                             # the data point to be explained\nxs = np.linspace(-2.0, 2.0, 201)\nwts = np.exp(-((xs - x0) ** 2) / (2 * 0.25 ** 2))  # proximity weights at x0\nA = np.c_[np.ones_like(xs), xs]                    # linear design (1, x)\nsw = np.sqrt(wts)\ncoef, *_ = np.linalg.lstsq(A * sw[:, None], h_hat(xs) * sw, rcond=None)\ng = A @ coef                                       # local linear surrogate\nnear = np.abs(xs - x0) < 0.3                       # around x0 ...\nfar = np.abs(xs + 1.5) < 0.3                       # ... vs far away\nerr_near = np.max(np.abs(h_hat(xs) - g)[near])\nerr_far = np.max(np.abs(h_hat(xs) - g)[far])\ncheck(\"[B-lime]  weighted linear fit matches h near x0, not globally\",\n      err_near < 0.3 and err_far > 4.0 * err_near)\n\n# few features: in d = 5, a hypothesis that depends on features 0 and 1\n# only yields a local surrogate with near-zero remaining coefficients\nZ = rng_l.standard_normal((400, 5)) * 0.3          # perturbations around a point\nz0 = np.array([0.8, -0.2, 0.5, 0.1, -0.4])\nZp = z0 + Z\nhz = np.sin(3.0 * Zp[:, 0]) + 0.5 * Zp[:, 1] ** 2  # depends on features 0, 1\nwz = np.exp(-np.sum(Z ** 2, axis=1) / (2 * 0.3 ** 2))\nAz = np.c_[np.ones(len(Zp)), Zp]\nswz = np.sqrt(wz)\ncz, *_ = np.linalg.lstsq(Az * swz[:, None], hz * swz, rcond=None)\nused = np.abs(cz[1:]) > 0.1                        # features the surrogate uses\ncheck(\"[B-lime]  the 5-feature local surrogate depends on 2 features\",\n      used[0] and used[1] and not used[2:].any())\n\n# write the LIME figure data (surrogate only inside the window)\ngwin = np.where(np.abs(xs - x0) <= 0.6, g, np.nan)\nwith open(OUT_DIR / \"interpretableml_lime.csv\", \"w\") as fh:\n    fh.write(\"x,h,g\\n\")\n    for xi, hi, gi in zip(xs, h_hat(xs), gwin):\n        fh.write(f\"{xi:.4f},{hi:.4f},{gi:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-eerm]** EERM idea (Zhang et al., 2024): the penalty charges the departure of the learned hypothesis from the predictions one user supplies. Those predictions are not elicited one data point at a time here; they are the predictions of a simpler proxy model the user considers interpretable, a linear-regression line. Adding lambda * mean((h(x) - userpred(x))^2) to the training of a degree-10 polynomial pulls the learned hypothesis toward that line; the EERM fit has (i) smaller departure from the user predictions and (ii) smaller test error than the fit without the penalty, at a slightly larger error on the trainset."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# EERM (Zhang et al., 2024): regularize the training of a high-capacity\n# model by penalizing the departure from the predictions one user supplies.\n# The user is characterized by those predictions; here they come from a\n# proxy model the user considers interpretable, rather than being elicited\n# one data point at a time.\nrng_e2 = np.random.default_rng(9)\nm2 = 14\nxt = np.sort(rng_e2.uniform(-1.0, 1.0, m2))        # training inputs\nf_true = lambda x: 1.2 * x + 0.4 * np.sin(2.0 * np.pi * x)  # noqa: E731\nyt = f_true(xt) + 0.4 * rng_e2.standard_normal(m2)\n\ndeg = 10                                           # high-capacity polynomial model\nPhi = np.vander(xt, deg + 1, increasing=True)      # train design\nxg = np.linspace(-1.0, 1.0, 201)                   # unlabeled grid for the penalty\nPsi = np.vander(xg, deg + 1, increasing=True)      # grid design\n\n# the user predictions, taken from a proxy model the user can follow:\n# the least-squares line through the data\ncl, *_ = np.linalg.lstsq(np.c_[np.ones(m2), xt], yt, rcond=None)\nsurr = cl[0] + cl[1] * xg                          # user predictions on the grid\n\nw_plain, *_ = np.linalg.lstsq(Phi, yt, rcond=None)  # unregularized training\nlam2 = 1.0                                          # EERM penalty weight\n# EERM objective (1/m)||y - Phi w||^2 + lam (1/N)||Psi w - surr||^2 as one\n# stacked least-squares problem\nAstack = np.vstack([Phi / np.sqrt(m2), np.sqrt(lam2 / len(xg)) * Psi])\nbstack = np.concatenate([yt / np.sqrt(m2), np.sqrt(lam2 / len(xg)) * surr])\nw_eerm, *_ = np.linalg.lstsq(Astack, bstack, rcond=None)\n\nh_plain, h_eerm = Psi @ w_plain, Psi @ w_eerm\ndisc_plain = np.mean((h_plain - surr) ** 2)        # departure from the user\ndisc_eerm = np.mean((h_eerm - surr) ** 2)\ntest_plain = np.mean((h_plain - f_true(xg)) ** 2)  # error vs the true mean\ntest_eerm = np.mean((h_eerm - f_true(xg)) ** 2)\ntr_plain = np.mean((yt - Phi @ w_plain) ** 2)\ntr_eerm = np.mean((yt - Phi @ w_eerm) ** 2)\ncheck(\"[B-eerm]  the penalty shrinks the departure from the user \"\n      \"predictions and the test error\",\n      disc_eerm < disc_plain and test_eerm < test_plain\n      and tr_plain <= tr_eerm)\n\nwith open(OUT_DIR / \"interpretableml_eerm.csv\", \"w\") as fh:\n    fh.write(\"x,hopaque,hsurr,heerm\\n\")\n    for r in zip(xg, h_plain, surr, h_eerm):\n        fh.write(\",\".join(f\"{v:.4f}\" for v in r) + \"\\n\")\nwith open(OUT_DIR / \"interpretableml_eermdata.csv\", \"w\") as fh:\n    fh.write(\"x,y\\n\")\n    for xi, yi in zip(xt, yt):\n        fh.write(f\"{xi:.4f},{yi:.4f}\\n\")\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass \"\n      f\"(train MSE dense {train_dense:.3f} < sparse {train_sparse:.3f}; \"\n      f\"test MSE sparse {test_sparse:.3f} < dense {test_dense:.3f}; \"\n      f\"EERM: disc {disc_plain:.3f} -> {disc_eerm:.3f}, \"\n      f\"test {test_plain:.3f} -> {test_eerm:.3f})\")\nprint(f\"wrote {OUT_DIR / 'interpretableml_lime.csv'}, \"\n      f\"{OUT_DIR / 'interpretableml_eerm.csv'}, \"\n      f\"{OUT_DIR / 'interpretableml_eermdata.csv'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}