{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "mean.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# mean \u2014 Python demo\n\nNumerical companion to the entry [mean](https://dictionaryofml.org/terms/mean.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nOne block per paragraph of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts. 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/mean.py`](https://dictionaryofml.org/terms/mean.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(), \"mean.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nmean.py \u2014 numerical companion to the glossary entry 'mean'.\n\nOne block per paragraph of the entry (marked [P...]): each block verifies\nnumerically what the corresponding statement asserts. Self-contained\n(numpy/matplotlib only), fixed seed.\n\nBlocks\n------\n[P-def]     The mean of a random vector is its expectation, the Lebesgue\n            integral of x with respect to its probability distribution:\n            the Monte Carlo average of iid draws converges to the exact\n            (analytic) mean as the sample grows.\n[P-sample]  A dataset D defines a discrete RV x~ = x^(I) with I uniform on\n            {1,...,m}; the mean of x~ equals the sample mean\n            (1/m) sum_r x^(r) exactly.\n[P-argmin]  For an RV with finite second moment, E{x} minimizes the risk\n            E{||x - c||^2}: the analytic mean beats every candidate c on a\n            grid, and the average loss over the sample has its minimum at\n            the sample mean.\n[P-erm]     Featureless regression: ERM with squared error loss over\n            labels y^(1..m) is solved by the sample mean \u2014 the closed-form\n            minimizer of (1/m) sum_r (y^(r) - h)^2 equals np.mean(y), and\n            every other h has larger average loss (setting of Fig. 1,\n            m = 5 labels with mean 4).\n[P-robust]  Outliers and the mean: moving one of m data points by delta\n            shifts the sample mean by exactly delta/m, without bound. With\n            at most gamma*m points corrupted, the gamma-trimmed mean and the\n            median stay bounded however far those points are moved.\n\nOutputs\n-------\nmean.png : preview figure (checking only).\n\nData generated by pythondemos/mean.py.\n\"\"\"\n\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\n\nrng = np.random.default_rng(42)\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-def]** The mean of a random vector is its expectation, the Lebesgue integral of x with respect to its probability distribution: the Monte Carlo average of iid draws converges to the exact (analytic) mean as the sample grows."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Mean = expectation = integral of x dP(x). Monte Carlo average of iid\n# draws from N(mu, C) converges to the analytic mean mu.\nprint(\"[P-def] Monte Carlo average converges to the analytic mean\")\nmu = np.array([1.0, -2.0])\nA = np.array([[1.0, 0.0], [0.5, 0.8]])\nerrs = []\nfor m in [10**2, 10**4, 10**6]:\n    x = rng.standard_normal((m, 2)) @ A.T + mu\n    errs.append(np.linalg.norm(x.mean(axis=0) - mu))\nprint(f\"    |MC mean - mu| for m=1e2,1e4,1e6: \"\n      f\"{errs[0]:.4f}, {errs[1]:.4f}, {errs[2]:.4f}\")\ncheck(\"MC error shrinks monotonically\", errs[0] > errs[1] > errs[2])\ncheck(\"MC mean at m=1e6 within 3e-3 of mu\", errs[2] < 3e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-sample]** A dataset D defines a discrete RV x~ = x^(I) with I uniform on {1,...,m}; the mean of x~ equals the sample mean (1/m) sum_r x^(r) exactly."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Discrete RV x~ = x^(I), I uniform on {1,...,m}: its mean is exactly the\n# sample mean of the dataset.\nprint(\"[P-sample] dataset-induced discrete RV has mean = sample mean\")\nD = rng.normal(size=(7, 2))\nprobs = np.full(7, 1 / 7)                      # P(I = r) = 1/m\nmean_discrete = (probs[:, None] * D).sum(axis=0)\ncheck(\"E{x^(I)} equals (1/m) sum_r x^(r)\",\n      np.allclose(mean_discrete, D.mean(axis=0)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-argmin]** For an RV with finite second moment, E{x} minimizes the risk E{||x - c||^2}: the analytic mean beats every candidate c on a grid, and the average loss over the sample has its minimum at the sample mean."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# E{x} = argmin_c E{||x - c||^2}. Empirically: risk(c) >= risk(mean) for\n# every c on a grid, with equality only at c = mean.\nprint(\"[P-argmin] the mean minimizes the expected squared distance\")\nx = rng.standard_normal((200000, 2)) @ A.T + mu\ndef emp_risk(c):\n    return np.mean(np.sum((x - c) ** 2, axis=1))\nrisk_mean = emp_risk(x.mean(axis=0))\ngrid = [x.mean(axis=0) + d for d in\n        (np.array([0.5, 0]), np.array([-0.3, 0.4]), np.array([0, -1.0]))]\ncheck(\"risk(mean) < risk(c) for all offset candidates c\",\n      all(emp_risk(c) > risk_mean for c in grid))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-erm]** Featureless regression: ERM with squared error loss over labels y^(1..m) is solved by the sample mean \u2014 the closed-form minimizer of (1/m) sum_r (y^(r) - h)^2 equals np.mean(y), and every other h has larger average loss (setting of Fig. 1, m = 5 labels with mean 4)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Featureless regression (Fig. 1 setting): labels 2,5,3,6,4; ERM with the\n# squared error loss is solved by the sample mean h_hat = 4.\nprint(\"[P-erm] featureless ERM with squared loss = sample mean\")\ny = np.array([2.0, 5.0, 3.0, 6.0, 4.0])\nh_grid = np.linspace(0, 8, 1601)\nemp = np.array([np.mean((y - h) ** 2) for h in h_grid])\nh_hat = h_grid[np.argmin(emp)]\ncheck(\"grid minimizer equals np.mean(y) = 4\", abs(h_hat - y.mean()) < 5e-3)\ncheck(\"sample mean is 4 (entry Fig. 1)\", np.isclose(y.mean(), 4.0))\ncheck(\"every other h has larger average loss\",\n      np.all(emp >= np.mean((y - y.mean()) ** 2) - 1e-12))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-robust]** Outliers and the mean: moving one of m data points by delta shifts the sample mean by exactly delta/m, without bound. With at most gamma*m points corrupted, the gamma-trimmed mean and the median stay bounded however far those points are moved."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# The entry's claim: moving one data point by delta shifts the sample mean\n# by delta/m, so a single corrupted point carries it arbitrarily far, while\n# the trimmed mean and the median do not move.\nprint(\"[P-robust] outliers move the mean without bound, the robust summaries stay bounded\")\n\nbase = rng.normal(size=101)\nm_pts = base.size\nfor delta in (10.0, 1e3, 1e6):\n    moved = base.copy()\n    moved[0] += delta\n    check(f\"one point moved by {delta:g}: mean shifts by delta/m\",\n          np.isclose(moved.mean() - base.mean(), delta / m_pts))\n\ndef trimmed(a, gamma):\n    \"\"\"Discard the floor(gamma*len) smallest and largest, average the rest.\"\"\"\n    k = int(np.floor(gamma * a.size))\n    return np.sort(a)[k: a.size - k].mean() if a.size - 2 * k > 0 else np.nan\n\nGAMMA = 0.1\nn_bad = int(np.floor(GAMMA * m_pts))          # at most gamma*m corrupted\n\ndef summaries(delta):\n    bad = base.copy()\n    bad[:n_bad] += delta\n    return bad.mean(), trimmed(bad, GAMMA), np.median(bad)\n\n# The robust summaries are NOT unchanged -- corrupting points changes which\n# order statistics are averaged -- but they stay BOUNDED as the corruption\n# grows without bound, which is the property that matters.\nrows = [summaries(d) for d in (1e3, 1e6, 1e9)]\nmean_shift = [abs(r[0] - base.mean()) for r in rows]\ntrim_shift = [abs(r[1] - trimmed(base, GAMMA)) for r in rows]\nmed_shift = [abs(r[2] - np.median(base)) for r in rows]\nprint(f\"    delta=1e3,1e6,1e9 -> mean moves by \"\n      f\"{mean_shift[0]:.3g}, {mean_shift[1]:.3g}, {mean_shift[2]:.3g}\")\nprint(f\"                        trimmed moves by \"\n      f\"{trim_shift[0]:.3g}, {trim_shift[1]:.3g}, {trim_shift[2]:.3g}\")\ncheck(f\"mean diverges with delta ({n_bad} of {m_pts} corrupted)\",\n      mean_shift[2] > 1e6)\ncheck(\"trimmed mean stays bounded as delta grows\",\n      max(trim_shift) < 1.0 and trim_shift[2] == trim_shift[0])\ncheck(\"median stays bounded as delta grows\",\n      max(med_shift) < 1.0 and med_shift[2] == med_shift[0])\n\n# ------------------------------------------------------------ preview\nfig, ax = plt.subplots(1, 2, figsize=(9, 3.2))\nax[0].loglog([1e2, 1e4, 1e6], errs, \"o-\")\nax[0].set_xlabel(\"m\"); ax[0].set_ylabel(\"|MC mean - mu|\")\nax[0].set_title(\"[P-def] average approaches the mean\")\nax[1].plot(h_grid, emp)\nax[1].axvline(y.mean(), ls=\"--\", c=\"k\")\nax[1].set_xlabel(\"h\"); ax[1].set_ylabel(\"average loss\")\nax[1].set_title(\"[P-erm] minimum at sample mean\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"mean.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}