{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "bootstrap.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# bootstrap \u2014 Python demo\n\nNumerical companion to the entry [bootstrap](https://dictionaryofml.org/terms/bootstrap.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nThree years of daily minimum and maximum temperature at Krems an der Donau stand in for the dataset D. The blocks below read those days as realizations of i.i.d. RVs, replace their unknown probability distribution by the empirical distribution of D, and draw from it. Self-contained (numpy/matplotlib only, data fetched from the GeoSphere Austria archive), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/bootstrap.py`](https://dictionaryofml.org/terms/bootstrap.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(), \"bootstrap.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nbootstrap.py \u2014 numerical companion to the glossary entry 'bootstrap'.\n\nThree years of daily minimum and maximum temperature at Krems an der\nDonau stand in for the dataset D. The blocks below read those days as\nrealizations of i.i.d. RVs, replace their unknown probability\ndistribution by the empirical distribution of D, and draw from it.\nSelf-contained (numpy/matplotlib only, data fetched from the GeoSphere\nAustria archive), fixed seed.\n\nBlocks\n------\n[B-data]    1096 days of (tmin, tmax) at Krems, 2022-2024. This cloud of\n            points is the dataset D; the entry's picture of it is the\n            scatterplot of the left panel.\n[B-empdist] The empirical distribution P^(D) puts mass 1/m on each of\n            the m days. Drawing m days from it is drawing from D with\n            replacement, and a draw holds about 0.63 m distinct days.\n            P^(D) is a complete distribution, so it can be sampled as\n            often as wanted: drawing 100 m days from it is no problem,\n            and the mean of such a draw converges to the mean of D, not\n            to the mean of the unknown P. The gain is replicates, not\n            information.\n[B-fit]     ERM with the squared error loss fits tmax = w0 + w1 tmin\n            on D, which is least squares. Refitting on\n            B = 500 bootstrap datasets gives B slopes whose 2.5th to\n            97.5th percentile range is a confidence interval for the\n            slope; it shrinks like 1/sqrt(m) when m is cut to a quarter.\n[B-testci]  A threshold rule predicts a frost night (tmin < 0) from the\n            day's tmax. The rule is learned once on a training set and\n            its accuracy measured on a held-out test set; resampling\n            that test set B times, with the rule held fixed, turns the\n            single accuracy into a confidence interval.\n[B-smooth]  Smearing each day with a Gaussian kernel turns the m spikes\n            of P^(D) into a density estimate, drawn as the contours of\n            the right panel. Sampling that density instead of the spikes\n            is the smoothed bootstrap: it produces days that D does not\n            contain, which the plain bootstrap never does.\n\nOutputs\n-------\nbootstrap_krems.csv    : tmin, tmax of the 1096 days (figure data).\nbootstrap_fit.csv      : the learned line as one two-point segment.\nbootstrap_lines.csv    : 25 bootstrap lines, two-point segments separated\n                         by empty lines.\nbootstrap_density.csv  : contour polylines of the smeared distribution,\n                         same separation.\nbootstrap.png          : preview figure (checking only).\n\nData generated by pythondemos/bootstrap.py.\n\"\"\"\n\nimport json\nimport urllib.request\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\nARCHIVE = \"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\nSTATION = 3805                      # Krems an der Donau\nNR_BOOTSTRAP = 500\nSEED = 20220101\n\nrng = np.random.default_rng(SEED)\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\ndef fetch(parameters, start, end):\n    \"\"\"Daily parameters of one station from the GeoSphere Austria archive.\"\"\"\n    url = (f\"{ARCHIVE}klima-v2-1d?parameters={','.join(parameters)}\"\n           f\"&station_ids={STATION}&start={start}&end={end}\")\n    with urllib.request.urlopen(url, timeout=180) as resp:\n        payload = json.load(resp)\n    values = payload[\"features\"][0][\"properties\"][\"parameters\"]\n    return np.stack([np.array(values[p][\"data\"], dtype=float)\n                     for p in parameters], axis=1)\n\n\ndef least_squares_line(x, y):\n    \"\"\"Slope and intercept of the least-squares fit of y on x.\"\"\"\n    design = np.stack([np.ones_like(x), x], axis=1)\n    weights = np.linalg.lstsq(design, y, rcond=None)[0]\n    return weights[1], weights[0]\n\n\ndef write_segments(name, header, segments):\n    \"\"\"Polylines, one blank line between them (pgfplots 'empty line=jump').\"\"\"\n    with open(OUT_DIR / name, \"w\") as f:\n        f.write(header + \"\\n\")\n        for k, seg in enumerate(segments):\n            if k:\n                f.write(\"\\n\")\n            for px, py in seg:\n                f.write(f\"{px:.4f},{py:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** 1096 days of (tmin, tmax) at Krems, 2022-2024. This cloud of points is the dataset D; the entry's picture of it is the scatterplot of the left panel."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "days = fetch([\"tlmin\", \"tlmax\"], \"2022-01-01\", \"2024-12-31\")\ntmin, tmax = days[:, 0], days[:, 1]\nm = tmin.size\nprint(f\"[B-data] {m} days at Krems an der Donau, 2022-2024; \"\n      f\"tmin {tmin.min():.1f} to {tmin.max():.1f}, \"\n      f\"tmax {tmax.min():.1f} to {tmax.max():.1f} degrees\")\ncheck(\"[B-data] three full years, no gaps\", m == 1096 and not np.isnan(days).any())\nwith open(OUT_DIR / \"bootstrap_krems.csv\", \"w\") as f:\n    f.write(\"tmin,tmax\\n\")\n    for a, b in zip(tmin, tmax):\n        f.write(f\"{a:.1f},{b:.1f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-empdist]** The empirical distribution P^(D) puts mass 1/m on each of the m days. Drawing m days from it is drawing from D with replacement, and a draw holds about 0.63 m distinct days. P^(D) is a complete distribution, so it can be sampled as often as wanted: drawing 100 m days from it is no problem, and the mean of such a draw converges to the mean of D, not to the mean of the unknown P. The gain is replicates, not information."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-empdist] P^(D) puts mass 1/m on each day; drawing from it is \"\n      \"drawing from D with replacement\")\nidx = rng.integers(0, m, size=m)\ncheck(\"[B-empdist] a draw of m days holds about 0.63 m distinct days\",\n      abs(np.unique(idx).size / m - (1 - np.exp(-1.0))) < 0.02)\ncheck(\"[B-empdist] every drawn day is a day of D\",\n      np.all(np.isin(tmin[idx], tmin)))\nhuge = rng.integers(0, m, size=100 * m)\ncheck(\"[B-empdist] a draw of 100 m days is no problem, and its mean \"\n      \"converges to the mean of D, not to anything new\",\n      abs(tmax[huge].mean() - tmax.mean()) < 0.1)\ncheck(\"[B-empdist] such a draw still contains no day outside D\",\n      np.unique(tmax[huge]).size <= np.unique(tmax).size)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fit]** ERM with the squared error loss fits tmax = w0 + w1 tmin on D, which is least squares. Refitting on B = 500 bootstrap datasets gives B slopes whose 2.5th to 97.5th percentile range is a confidence interval for the slope; it shrinks like 1/sqrt(m) when m is cut to a quarter."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "slope, intercept = least_squares_line(tmin, tmax)\nprint(f\"[B-fit] tmax = {intercept:.2f} + {slope:.2f} tmin on all {m} days\")\nslopes = np.empty(NR_BOOTSTRAP)\nlines = []\nfor b in range(NR_BOOTSTRAP):\n    draw = rng.integers(0, m, size=m)\n    slopes[b], icept = least_squares_line(tmin[draw], tmax[draw])\n    if b < 25:\n        lines.append([(tmin.min(), icept + slopes[b] * tmin.min()),\n                      (tmin.max(), icept + slopes[b] * tmin.max())])\nlo, hi = np.percentile(slopes, [2.5, 97.5])\nprint(f\"[B-fit] 95% confidence interval for the slope: [{lo:.3f}, {hi:.3f}]\")\ncheck(\"[B-fit] the interval covers the slope learned from all of D\",\n      lo <= slope <= hi)\ncheck(\"[B-fit] the interval is narrow, a few percent of the slope\",\n      (hi - lo) / slope < 0.12)\nquarter = rng.choice(m, size=m // 4, replace=False)\nnarrow = np.array([least_squares_line(*(lambda d: (tmin[quarter][d],\n                                                   tmax[quarter][d]))(\n    rng.integers(0, m // 4, size=m // 4)))[0] for _ in range(NR_BOOTSTRAP)])\nwide = np.percentile(narrow, 97.5) - np.percentile(narrow, 2.5)\nprint(f\"[B-fit] on a quarter of the days the interval is {wide / (hi - lo):.1f} \"\n      \"times as wide\")\ncheck(\"[B-fit] cutting m to a quarter roughly doubles the interval\",\n      1.6 < wide / (hi - lo) < 2.6)\n# the learned line goes in its own file: the book figure draws it thick and\n# the replicates thin, so one file with one style would hide the difference\nwrite_segments(\"bootstrap_fit.csv\", \"tmin,tmax\",\n               [[(tmin.min(), intercept + slope * tmin.min()),\n                 (tmin.max(), intercept + slope * tmin.max())]])\nwrite_segments(\"bootstrap_lines.csv\", \"tmin,tmax\", lines)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-testci]** A threshold rule predicts a frost night (tmin < 0) from the day's tmax. The rule is learned once on a training set and its accuracy measured on a held-out test set; resampling that test set B times, with the rule held fixed, turns the single accuracy into a confidence interval."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "frost = (tmin < 0.0).astype(int)\norder = rng.permutation(m)\ntrain, test = order[:800], order[800:]\ngrid = np.linspace(tmax.min(), tmax.max(), 400)\nerrors = [(np.mean((tmax[train] <= t).astype(int) != frost[train]), t)\n          for t in grid]\nthreshold = min(errors)[1]\npredict = lambda x: (x <= threshold).astype(int)\nacc_test = float(np.mean(predict(tmax[test]) == frost[test]))\nprint(f\"[B-testci] threshold {threshold:.1f} degrees learned on \"\n      f\"{train.size} days; accuracy on the {test.size} test days \"\n      f\"{acc_test:.3f}\")\ncheck(\"[B-testci] the learned rule is better than always predicting \"\n      \"no frost\", acc_test > 1 - frost[test].mean())\naccs = np.array([float(np.mean(predict(tmax[test][d]) == frost[test][d]))\n                 for d in rng.integers(0, test.size,\n                                       size=(NR_BOOTSTRAP, test.size))])\nacc_lo, acc_hi = np.percentile(accs, [2.5, 97.5])\nprint(f\"[B-testci] 95% confidence interval for that accuracy: \"\n      f\"[{acc_lo:.3f}, {acc_hi:.3f}]\")\ncheck(\"[B-testci] the interval covers the measured accuracy\",\n      acc_lo <= acc_test <= acc_hi)\ncheck(\"[B-testci] it is several points wide, which the single number \"\n      \"does not show\", 0.02 < acc_hi - acc_lo < 0.12)\ncheck(\"[B-testci] the hypothesis never changes: only the test set is \"\n      \"redrawn\", predict(np.array([threshold])).item() == 1)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-smooth]** Smearing each day with a Gaussian kernel turns the m spikes of P^(D) into a density estimate, drawn as the contours of the right panel. Sampling that density instead of the spikes is the smoothed bootstrap: it produces days that D does not contain, which the plain bootstrap never does."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "h_min, h_max = 1.4, 1.8                       # kernel widths, degrees\ngx = np.linspace(tmin.min() - 3, tmin.max() + 3, 160)\ngy = np.linspace(tmax.min() - 3, tmax.max() + 3, 160)\nGX, GY = np.meshgrid(gx, gy)\ndens = np.zeros_like(GX)\nfor a, b in zip(tmin, tmax):\n    dens += np.exp(-0.5 * (((GX - a) / h_min) ** 2 + ((GY - b) / h_max) ** 2))\ndens /= m * 2 * np.pi * h_min * h_max\ncell = (gx[1] - gx[0]) * (gy[1] - gy[0])\nprint(f\"[B-smooth] kernel widths {h_min} and {h_max} degrees; the smeared \"\n      f\"density integrates to {dens.sum() * cell:.3f}\")\ncheck(\"[B-smooth] the density integrates to one\", abs(dens.sum() * cell - 1) < 0.02)\ncounts, edges = np.histogram(tmin, bins=20)\nmode_tmin = 0.5 * (edges[counts.argmax()] + edges[counts.argmax() + 1])\npeak_tmin = gx[np.unravel_index(dens.argmax(), dens.shape)[1]]\ncheck(f\"[B-smooth] its peak ({peak_tmin:.1f} degrees) sits where the days \"\n      f\"pile up ({mode_tmin:.1f}), not at the median ({np.median(tmin):.1f})\",\n      abs(peak_tmin - mode_tmin) < 1.5)\nsmoothed = (tmin[rng.integers(0, m, size=20000)]\n            + h_min * rng.standard_normal(20000))\ncheck(\"[B-smooth] the smoothed bootstrap produces days outside D, which \"\n      \"the plain bootstrap never does\",\n      not np.any(np.isin(smoothed, tmin)))\nlevels = np.array([0.1, 0.4, 0.7]) * dens.max()\n\n# ------------------------------------------------------------- preview\nfig, ax = plt.subplots(1, 3, figsize=(13.2, 3.6))\nax[0].plot(tmin, tmax, \".\", color=\"0.55\", markersize=2.5, label=\"day\")\nfor seg in lines[:25]:\n    ax[0].plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]],\n               \"-\", color=\"0.3\", linewidth=0.4)\nax[0].plot([tmin.min(), tmin.max()],\n           [intercept + slope * tmin.min(), intercept + slope * tmin.max()],\n           \"k-\", linewidth=1.8, label=\"learned line\")\nax[0].plot([], [], \"-\", color=\"0.3\", linewidth=0.4, label=\"bootstrap lines\")\nax[0].set_xlabel(\"minimum temperature (degC)\")\nax[0].set_ylabel(\"maximum temperature (degC)\")\nax[0].set_title(\"[B-fit] the dataset and the spread of the fits\", fontsize=10)\nax[0].legend(frameon=False, fontsize=8, loc=\"upper left\")\nax[1].plot(tmin, tmax, \".\", color=\"0.75\", markersize=2)\ncontours = ax[1].contour(GX, GY, dens, levels=levels, colors=\"black\",\n                         linewidths=0.9)\n# the book figure draws the same contours, so its polylines are taken from\n# this panel rather than from a second, throwaway axes\nsegments = [seg for level in contours.allsegs for seg in level if len(seg) > 8]\nwrite_segments(\"bootstrap_density.csv\", \"tmin,tmax\", segments)\nprint(f\"[B-smooth] {len(segments)} contour polylines written\")\nax[1].set_xlabel(\"minimum temperature (degC)\")\nax[1].set_ylabel(\"maximum temperature (degC)\")\nax[1].set_title(\"[B-smooth] the same days, smeared into a density\",\n                fontsize=10)\nax[2].hist(slopes, bins=30, color=\"white\", edgecolor=\"black\")\nax[2].axvline(slope, color=\"black\", linestyle=\"-\", linewidth=1.5,\n              label=\"slope on all days\")\nax[2].axvline(lo, color=\"black\", linestyle=\"--\", linewidth=1.0,\n              label=\"95% interval\")\nax[2].axvline(hi, color=\"black\", linestyle=\"--\", linewidth=1.0)\nax[2].set_xlabel(\"slope of the bootstrap line\")\nax[2].set_ylabel(\"number of bootstrap datasets\")\nax[2].set_title(\"[B-fit] the slopes of 500 replicates\", fontsize=10)\nax[2].legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"bootstrap.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}