{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "randomforest.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# random forest \u2014 Python demo\n\nNumerical companion to the entry [random forest](https://dictionaryofml.org/terms/randomforest.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nNumerical companion to the glossary entry 'randomforest' (random forest). The GeoSphere Austria weather station Krems (station id 3805) records the minimum and maximum air temperature of each day; this script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them to randomforest_weather.csv. Each day is a data point whose feature is the minimum temperature tmin and whose label is the maximum temperature tmax.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/randomforest.py`](https://dictionaryofml.org/terms/randomforest.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(), \"randomforest.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Three bootstrap regression trees and their average at Krems: the\nforest fits no worse than its trees on average.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'randomforest' (random\nforest).  The GeoSphere Austria weather station Krems (station id 3805)\nrecords the minimum and maximum air temperature of each day; this\nscript downloads the records for 2024 from the GeoSphere data hub\n(dataset klima-v2-1d) and writes them to randomforest_weather.csv.\nEach day is a data point whose feature is the minimum temperature tmin\nand whose label is the maximum temperature tmax.\n\nThree depth-2 regression trees are trained, each on its own bootstrap\nsample of the 366 days, and the random forest predicts by averaging\nthe three tree predictions.  With a single feature there is no room\nfor feature subsampling, so the forest's randomness comes from the\nbootstrap alone, as in bagging.\n\nThe demo checks the entry's central claims: each tree is a piecewise\nconstant map with few pieces; the bootstrap makes the trees differ;\nand the squared-error risk of the averaged prediction on the full\ndataset never exceeds the average of the trees' squared-error risks\n(the algebraic identity behind variance reduction by averaging).\n\nDeterministic: the bootstrap samples use fixed seeds.  Self-contained:\nnumpy + matplotlib only (stdlib urllib for the download).\n\nBlocks\n------\n[B-fetch] Download the 366 daily temperature pairs at Krems for 2024\n          and write them to randomforest_weather.csv; check the count\n          and one pinned value against the archive.\n[B-trees] Train three depth-2 regression trees, each on a bootstrap\n          sample (fixed seeds); check that each tree is piecewise\n          constant with at most four pieces, fits its bootstrap sample\n          better than a constant, and that the trees differ pairwise.\n[B-forest] Average the three trees into the random forest prediction;\n          check that the forest curve is the pointwise mean of the\n          tree curves and that the forest's squared-error risk on all\n          366 days is smaller than the average of the trees' risks.\n\nOutputs\n-------\nrandomforest_weather.csv : date, tmin, tmax for the 366 days of 2024\nrandomforest_tree1.csv   : tmin, pred -- first tree on a tmin grid\nrandomforest_tree2.csv   : tmin, pred -- second tree on the grid\nrandomforest_tree3.csv   : tmin, pred -- third tree on the grid\nrandomforest_forest.csv  : tmin, pred -- the averaged (forest) curve\nrandomforest.png         : preview (checking only) -- scatterplot of\n                           the days, the three tree curves, and the\n                           forest curve\n\"\"\"\n\nimport json\nimport urllib.request\nfrom pathlib import Path\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\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}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fetch]** Download the 366 daily temperature pairs at Krems for 2024 and write them to randomforest_weather.csv; check the count and one pinned value against the archive."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "URL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       \"klima-v2-1d?parameters=tlmin,tlmax&station_ids=3805\"\n       \"&start=2024-01-01&end=2024-12-31\")\nwith urllib.request.urlopen(URL, timeout=120) as resp:\n    payload = json.load(resp)\nparams = payload[\"features\"][0][\"properties\"][\"parameters\"]\nstamps = [t[:10] for t in payload[\"timestamps\"]]\ntmin = np.array(params[\"tlmin\"][\"data\"], dtype=float)\ntmax = np.array(params[\"tlmax\"][\"data\"], dtype=float)\nwith open(OUT_DIR / \"randomforest_weather.csv\", \"w\") as f:\n    f.write(\"date,tmin,tmax\\n\")\n    for d, a, b in zip(stamps, tmin, tmax):\n        f.write(f\"{d},{a:.1f},{b:.1f}\\n\")\ncheck(\"[B-fetch] 366 daily temperature pairs downloaded for 2024\",\n      len(tmin) == 366)\nfeb1 = stamps.index(\"2024-02-01\")\ncheck(\"[B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4)\",\n      tmin[feb1] == -3.8 and tmax[feb1] == 10.4)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-trees]** Train three depth-2 regression trees, each on a bootstrap sample (fixed seeds); check that each tree is piecewise constant with at most four pieces, fits its bootstrap sample better than a constant, and that the trees differ pairwise."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def best_split(x, y):\n    \"\"\"Threshold minimizing the summed squared error of the two parts,\n    or None if no split with at least 5 data points per part exists.\"\"\"\n    order = np.argsort(x)\n    xs, ys = x[order], y[order]\n    best, best_sse = None, np.inf\n    for i in range(5, len(xs) - 5 + 1):\n        if xs[i - 1] == xs[i]:\n            continue\n        left, right = ys[:i], ys[i:]\n        sse = ((left - left.mean()) ** 2).sum() \\\n            + ((right - right.mean()) ** 2).sum()\n        if sse < best_sse:\n            best_sse, best = sse, (xs[i - 1] + xs[i]) / 2\n    return best\n\n\ndef depth2_tree(x, y):\n    \"\"\"Thresholds and leaf means of a depth-2 regression tree.\"\"\"\n    t0 = best_split(x, y)\n    cuts = []\n    for side in (x <= t0, x > t0):\n        t = best_split(x[side], y[side])\n        if t is not None:\n            cuts.append(t)\n    cuts = sorted(cuts + [t0])\n    edges = [-np.inf] + cuts + [np.inf]\n    means = [y[(x > lo) & (x <= hi)].mean()\n             for lo, hi in zip(edges[:-1], edges[1:])]\n    return np.array(cuts), np.array(means)\n\n\ndef predict(cuts, means, xq):\n    return means[np.searchsorted(cuts, xq)]\n\n\ngrid = np.linspace(tmin.min() - 1.0, tmin.max() + 1.0, 400)\ntrees, curves = [], []\nfor seed in (1, 2, 3):\n    rng = np.random.default_rng(seed)\n    idx = rng.integers(0, len(tmin), len(tmin))     # bootstrap sample\n    cuts, means = depth2_tree(tmin[idx], tmax[idx])\n    trees.append((idx, cuts, means))\n    curves.append(predict(cuts, means, grid))\n    np.savetxt(OUT_DIR / f\"randomforest_tree{seed}.csv\",\n               np.stack([grid, curves[-1]], 1),\n               delimiter=\",\", header=\"tmin,pred\", comments=\"\", fmt=\"%.2f\")\ncheck(\"[B-trees] each tree is piecewise constant with at most 4 pieces\",\n      all(len(means) <= 4 for _, _, means in trees))\ncheck(\"[B-trees] each tree fits its bootstrap sample better than a constant\",\n      all(((tmax[idx] - predict(cuts, means, tmin[idx])) ** 2).mean()\n          < tmax[idx].var() for idx, cuts, means in trees))\ncheck(\"[B-trees] the bootstrap makes the trees differ pairwise\",\n      all(np.max(np.abs(curves[i] - curves[j])) > 0.5\n          for i in range(3) for j in range(i + 1, 3)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-forest]** Average the three trees into the random forest prediction; check that the forest curve is the pointwise mean of the tree curves and that the forest's squared-error risk on all 366 days is smaller than the average of the trees' risks."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "forest = np.mean(curves, axis=0)\nnp.savetxt(OUT_DIR / \"randomforest_forest.csv\",\n           np.stack([grid, forest], 1),\n           delimiter=\",\", header=\"tmin,pred\", comments=\"\", fmt=\"%.2f\")\ntree_risks = [((tmax - predict(cuts, means, tmin)) ** 2).mean()\n              for _, cuts, means in trees]\nforest_pred = np.mean([predict(cuts, means, tmin)\n                       for _, cuts, means in trees], axis=0)\nforest_risk = ((tmax - forest_pred) ** 2).mean()\ncheck(\"[B-forest] the forest curve is the pointwise mean of the tree curves\",\n      np.allclose(forest, np.mean(curves, axis=0)))\ncheck(\"[B-forest] the averaged prediction has smaller squared-error risk \"\n      \"than the trees on average\",\n      forest_risk < np.mean(tree_risks) - 1e-9)\nprint(f\"  tree risks {[round(r, 2) for r in tree_risks]}, \"\n      f\"forest risk {forest_risk:.2f}\")\n\nfig, ax = plt.subplots(figsize=(6.4, 4.4))\nax.plot(tmin, tmax, \"o\", ms=2, color=\"0.6\", label=\"days of 2024 (Krems)\")\nfor c, (style, name) in zip(curves, [(\":\", \"tree 1\"), (\"--\", \"tree 2\"),\n                                     (\"-.\", \"tree 3\")]):\n    ax.plot(grid, c, style, lw=1.2, label=name)\nax.plot(grid, forest, \"-\", lw=2.2, color=\"black\", label=\"random forest\")\nax.set_xlabel(\"minimum temperature of the day (\u00b0C)\")\nax.set_ylabel(\"maximum temperature of the day (\u00b0C)\")\nax.set_title(\"three bootstrap trees and their average (random forest)\")\nax.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"randomforest.png\", dpi=150)\n\nfailed = [name for name, ok in report if not ok]\nprint(f\"{len(report) - len(failed)}/{len(report)} checks passed\"\n      + (f\", FAILED: {failed}\" if failed else \"\"))"
  }
 ]
}