{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "ensemble.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# ensemble \u2014 Python demo\n\nNumerical companion to the entry [ensemble](https://dictionaryofml.org/terms/ensemble.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 quantitative claims. First, a majority vote of three different classifiers over patches of an aerial photograph (vineyard or not) misclassifies fewer patches than any single one of them. Second, the average of B identically distributed predictions with variance sigma^2 and pairwise correlation rho has variance rho sigma^2 + (1 - rho) sigma^2 / B: the second term vanishes as B grows, the first remains. Self-contained (numpy/matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/ensemble.py`](https://dictionaryofml.org/terms/ensemble.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(), \"ensemble.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nensemble.py \u2014 numerical companion to the glossary entry 'ensemble'.\n\nPurpose\n-------\nBacks the entry's two quantitative claims.  First, a majority vote of\nthree different classifiers over patches of an aerial photograph\n(vineyard or not) misclassifies fewer patches than any single one of\nthem.  Second, the average of B identically distributed predictions with\nvariance sigma^2 and pairwise correlation rho has variance\nrho sigma^2 + (1 - rho) sigma^2 / B: the second term vanishes as B grows,\nthe first remains.  Self-contained (numpy/matplotlib only), fixed seed.\n\nSetup\n-----\nPatches: 600 labelled patches with two features, a color score and a\ntexture score, drawn around two overlapping centers (vineyard / not);\nthe label is flipped for 8 percent of the patches.  300 patches form\nthe training set, the other 300 are kept for measuring the misclassified\nfraction.  Three base learners, each a decision tree with two levels of\nsplits, each trained on its own random subset of 120 training patches,\nso that their errors are only partly shared.  The aggregation rule is\nthe majority vote.\n\nBlocks\n------\n[B-vote]  Over 200 repetitions with fresh random subsets, the majority\n          vote misclassifies fewer of the 300 held-back patches, on\n          average, than a single base learner, and it beats the average\n          of its three base learners in at least 80 percent of the\n          repetitions.\n[B-var]   For sigma^2 = 1 and rho = 0.3, the empirical variance of the\n          average of B predictions, over 20000 draws, matches\n          rho sigma^2 + (1 - rho) sigma^2 / B within 5 percent for every\n          B in {1, 2, 3, 5, 10, 20, 50}.\n[B-limit] At B = 50 the variance of the average is within 10 percent of\n          rho sigma^2 = 0.3: the correlated part is what averaging\n          cannot remove.\n\nOutputs\n-------\nensemble_vote.csv     : one row per repetition: misclassified fraction\n                        of the three base learners and of the majority\n                        vote.\nensemble_variance.csv : B, empirical variance of the average, and the\n                        formula value.\nensemble.png          : matplotlib preview (checking only).\n\"\"\"\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\nrng = np.random.default_rng(0)\n\n# ---------------------------------------------------------------- patches\nn = 600\ny = rng.integers(0, 2, n)                                  # 1 = vineyard\ncenters = np.array([[0.35, 0.35], [0.65, 0.65]])\nX = centers[y] + 0.17 * rng.standard_normal((n, 2))       # color, texture\nflip = rng.random(n) < 0.08\ny = np.where(flip, 1 - y, y)\ntr = np.arange(n) < 300\nX_tr, y_tr, X_te, y_te = X[tr], y[tr], X[~tr], y[~tr]\n\n\ndef best_split(Xs, ys, j):\n    \"\"\"Threshold on feature j that misclassifies fewest patches, with the\n    class assigned to each side; returns (threshold, left class, right class,\n    misclassified count).\"\"\"\n    order = np.argsort(Xs[:, j]); v = Xs[order, j]; lab = ys[order]\n    best = (None, 0, 1, len(ys) + 1)\n    for k in range(1, len(v)):\n        if v[k] == v[k - 1]:\n            continue\n        thr = 0.5 * (v[k] + v[k - 1])\n        left, right = lab[:k], lab[k:]\n        lc = int(left.mean() >= 0.5); rc = int(right.mean() >= 0.5)\n        err = int((left != lc).sum() + (right != rc).sum())\n        if err < best[3]:\n            best = (thr, lc, rc, err)\n    return best\n\n\ndef stump(Xs, ys, j):\n    thr, lc, rc, _ = best_split(Xs, ys, j)\n    return lambda Z: np.where(Z[:, j] <= thr, lc, rc)\n\n\ndef tree2(Xs, ys):\n    \"\"\"Two levels of splits: the better feature at the root, then the\n    better feature on each side.\"\"\"\n    cands = [best_split(Xs, ys, j) for j in (0, 1)]\n    j0 = int(np.argmin([c[3] for c in cands])); thr0 = cands[j0][0]\n    left = Xs[:, j0] <= thr0\n    subs = []\n    for side in (left, ~left):\n        cs = [best_split(Xs[side], ys[side], j) for j in (0, 1)]\n        j1 = int(np.argmin([c[3] for c in cs])); subs.append((j1, cs[j1]))\n\n    def predict(Z):\n        out = np.empty(len(Z), dtype=int)\n        for side, (j1, (thr, lc, rc, _)) in zip((Z[:, j0] <= thr0, Z[:, j0] > thr0), subs):\n            out[side] = np.where(Z[side, j1] <= thr, lc, rc)\n        return out\n    return predict"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-vote]** Over 200 repetitions with fresh random subsets, the majority vote misclassifies fewer of the 300 held-back patches, on average, than a single base learner, and it beats the average of its three base learners in at least 80 percent of the repetitions."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "reps = 200\nrows = []\nfor _ in range(reps):\n    learners = []\n    for _ in range(3):\n        idx = rng.choice(300, 120, replace=False)          # own random subset\n        learners.append(tree2(X_tr[idx], y_tr[idx]))\n    votes = np.stack([h(X_te) for h in learners])          # 3 x 300\n    errs = [float((v != y_te).mean()) for v in votes]\n    vote = (votes.sum(axis=0) >= 2).astype(int)\n    rows.append(errs + [float((vote != y_te).mean())])\nrows = np.array(rows)\nmean_single, mean_vote = float(rows[:, :3].mean()), float(rows[:, 3].mean())\nwins = float((rows[:, 3] < rows[:, :3].mean(axis=1)).mean())\ncheck(f\"[B-vote]  misclassified fraction over {reps} repetitions: single base \"\n      f\"learner {mean_single:.3f} on average, majority vote {mean_vote:.3f}; the \"\n      f\"vote beats its base learners' average in {100 * wins:.0f}% of them\",\n      mean_vote < mean_single and wins >= 0.8)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-var]** For sigma^2 = 1 and rho = 0.3, the empirical variance of the average of B predictions, over 20000 draws, matches rho sigma^2 + (1 - rho) sigma^2 / B within 5 percent for every B in {1, 2, 3, 5, 10, 20, 50}."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "sigma2, rho, draws = 1.0, 0.3, 20000\nBs = [1, 2, 3, 5, 10, 20, 50]\nemp, formula = [], []\nfor B in Bs:\n    shared = rng.standard_normal((draws, 1)) * np.sqrt(rho * sigma2)\n    own = rng.standard_normal((draws, B)) * np.sqrt((1 - rho) * sigma2)\n    avg = (shared + own).mean(axis=1)                      # correlation rho\n    emp.append(float(avg.var())); formula.append(rho * sigma2 + (1 - rho) * sigma2 / B)\ndev = max(abs(e - f) / f for e, f in zip(emp, formula))\ncheck(f\"[B-var]   variance of the average vs rho s^2 + (1-rho) s^2/B: \"\n      f\"largest relative deviation {dev:.3f}\", dev < 0.05)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-limit]** At B = 50 the variance of the average is within 10 percent of rho sigma^2 = 0.3: the correlated part is what averaging cannot remove."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "check(f\"[B-limit] at B=50 the variance is {emp[-1]:.3f}, the correlated part \"\n      f\"rho s^2 = {rho * sigma2:.1f} remains\", abs(emp[-1] - rho * sigma2) < 0.1 * rho * sigma2)\n\n# ---------------------------------------------------------------- CSV\nwith open(OUT_DIR / \"ensemble_vote.csv\", \"w\") as fh:\n    fh.write(\"rep,tree1,tree2,tree3,vote\\n\")\n    for k, r in enumerate(rows):\n        fh.write(f\"{k},{r[0]:.4f},{r[1]:.4f},{r[2]:.4f},{r[3]:.4f}\\n\")\nwith open(OUT_DIR / \"ensemble_variance.csv\", \"w\") as fh:\n    fh.write(\"B,empirical,formula\\n\")\n    for B, e, f in zip(Bs, emp, formula):\n        fh.write(f\"{B},{e:.4f},{f:.4f}\\n\")\n\n# -------------------------------------------------------------- preview\nfig, (ax, ax2) = plt.subplots(1, 2, figsize=(9.4, 3.8))\nax.plot(rows[:, :3].mean(axis=1), rows[:, 3], \"k.\", ms=4, label=\"one repetition\")\nlim = [rows.min() - 0.01, rows.max() + 0.01]\nax.plot(lim, lim, \"k--\", lw=0.8, label=\"equal error\")\nax.set_xlabel(\"average misclassified fraction of the three base learners\")\nax.set_ylabel(\"misclassified fraction of the majority vote\")\nax.set_title(\"majority vote vs its base learners (200 repetitions)\")\nax.legend(frameon=False, fontsize=8)\nax2.plot(Bs, formula, \"k-\", lw=1.3, label=r\"$\\rho\\sigma^2 + (1-\\rho)\\sigma^2/B$\")\nax2.plot(Bs, emp, \"ko\", mfc=\"none\", ms=6, label=\"empirical variance of the average\")\nax2.axhline(rho * sigma2, color=\"k\", ls=\":\", lw=1, label=r\"$\\rho\\sigma^2$ (shared part)\")\nax2.set_xscale(\"log\"); ax2.set_xlabel(\"number of base learners B\"); ax2.set_ylabel(\"variance of the average prediction\")\nax2.set_title(r\"averaging keeps the shared part ($\\rho = 0.3$)\")\nax2.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"ensemble.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(f\"wrote {OUT_DIR / 'ensemble_vote.csv'}, {OUT_DIR / 'ensemble_variance.csv'}, \"\n      f\"{OUT_DIR / 'ensemble.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}