{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "probdist.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# probability distribution \u2014 Python demo\n\nNumerical companion to the entry [probability distribution](https://dictionaryofml.org/terms/probdist.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/probdist.py`](https://dictionaryofml.org/terms/probdist.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(), \"probdist.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nprobdist.py \u2014 numerical companion to the glossary entry\n'probability distribution'.\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\nThe method used throughout the first blocks is the simplest one\navailable: it reads a data set of numbers and delivers their average.\nThat keeps every quantity below available in closed form, so a check\ncompares a measured number against a formula rather than against\nanother simulation.\n\nBlocks\n------\n[P-method]  One distribution, many data sets, a different output each\n            time: the outputs of the method are realizations of another\n            RV, whose spread shrinks as sigma/sqrt(m) with the data set\n            size m. The distribution of that output RV depends on the\n            common distribution P and on the method: a different method\n            (taking the middle value of the sorted data set) run on the\n            same data sets has a larger output spread.\n[P-typical] The distribution decides which data points are typical\n            (relative frequencies of a growing sample converge to it),\n            and it can be estimated from a data set: the empirical\n            frequencies computed from one large data set recover the\n            distribution that generated it.\n[P-specify] How a distribution is specified: a binary RV by the single\n            probability P(y = 0), and a continuous real-valued RV by a\n            pdf p, for which P(x in [a, b]) ~ p(a)|b - a| on a short\n            interval.\n[P-picture] The distribution visualized: the joint pdf of a data point\n            (feature, label) as a grayscale over the plane, a\n            two-component Gaussian mixture. Paralleling the entry's\n            Fig. 1, three data sets of realizations of iid RVs with\n            this distribution are drawn, and three hypothesis maps are\n            learned from them by the same polynomial regression method.\n            The checks verify that the plotted window carries the\n            probability mass, that the data sets concentrate where the\n            pdf is large, that the three learned maps differ, and that\n            they differ most where the pdf is small.\n\nOutputs\n-------\nprobdist.png             : preview figure (checking only).\nprobdist_density.csv     : joint pdf on a grid (columns x, y, p; row-wise\n                           in y).\nprobdist_train1.csv, probdist_train2.csv, probdist_train3.csv\n                         : the three data sets (columns x, y).\nprobdist_hypotheses.csv  : the three learned hypothesis curves\n                           (columns x, yhat1, yhat2, yhat3).\n\nData generated by pythondemos/probdist.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\nMU, SIGMA = 1.0, 2.0            # the common distribution: N(MU, SIGMA^2)\nNR_DATASETS = 4000\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\ndef method(data):\n    \"\"\"The ML method: it delivers the average of the data set.\"\"\"\n    return float(np.mean(data))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-method]** One distribution, many data sets, a different output each time: the outputs of the method are realizations of another RV, whose spread shrinks as sigma/sqrt(m) with the data set size m. The distribution of that output RV depends on the common distribution P and on the method: a different method (taking the middle value of the sorted data set) run on the same data sets has a larger output spread."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-method] one distribution, many data sets, a different output each time\")\noutputs = {}\nfor m in (25, 400):\n    outputs[m] = np.array([method(rng.normal(MU, SIGMA, m))\n                           for _ in range(NR_DATASETS)])\nspread = {m: float(np.std(out)) for m, out in outputs.items()}\nprint(f\"    output spread: {spread[25]:.4f} at m=25, {spread[400]:.4f} at m=400\"\n      f\"  (sigma/sqrt(m) = {SIGMA / np.sqrt(25):.4f}, \"\n      f\"{SIGMA / np.sqrt(400):.4f})\")\ncheck(\"two data sets from the same distribution give different outputs\",\n      outputs[25][0] != outputs[25][1])\ncheck(\"the spread of the outputs shrinks with the data set size\",\n      spread[400] < spread[25])\nfor m in (25, 400):\n    check(f\"the spread matches sigma/sqrt(m) within 5% at m={m}\",\n          abs(spread[m] - SIGMA / np.sqrt(m)) / (SIGMA / np.sqrt(m)) < 0.05)\n# the distribution of the output RV depends on P and on the method: the\n# middle value of a sorted Gaussian data set has a larger spread than\n# the average\nmed_outputs = np.array([float(np.median(rng.normal(MU, SIGMA, 25)))\n                        for _ in range(NR_DATASETS)])\nprint(f\"    output spread at m=25: {spread[25]:.4f} (average), \"\n      f\"{float(np.std(med_outputs)):.4f} (middle value)\")\ncheck(\"a different method (the middle value) has a different output \"\n      \"distribution: its spread is larger\",\n      float(np.std(med_outputs)) > 1.1 * spread[25])"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-typical]** The distribution decides which data points are typical (relative frequencies of a growing sample converge to it), and it can be estimated from a data set: the empirical frequencies computed from one large data set recover the distribution that generated it."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-typical] typical data points; the distribution estimated from a data set\")\np_true = np.array([0.5, 0.3, 0.2])              # distribution on {0, 1, 2}\nerrs = []\nfor m in (10**2, 10**4, 10**6):\n    draws = rng.choice(3, size=m, p=p_true)\n    errs.append(np.max(np.abs(np.bincount(draws, minlength=3) / m - p_true)))\nprint(f\"    max |frequency - p| for m=1e2,1e4,1e6: \"\n      f\"{errs[0]:.4f}, {errs[1]:.4f}, {errs[2]:.4f}\")\ncheck(\"relative frequencies converge to the distribution\", errs[0] > errs[2])\ncheck(\"the empirical frequencies of one large data set estimate the \"\n      \"distribution that generated it (within 2e-3)\", errs[2] < 2e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-specify]** How a distribution is specified: a binary RV by the single probability P(y = 0), and a continuous real-valued RV by a pdf p, for which P(x in [a, b]) ~ p(a)|b - a| on a short interval."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-specify] one probability specifies a binary RV; a pdf a continuous one\")\np0 = 0.73\ny = (rng.uniform(size=10**6) >= p0).astype(int)   # P(y = 0) = p0\nf0 = float(np.mean(y == 0))\ncheck(\"empirical P(y = 0) recovers p0 = 0.73\", abs(f0 - p0) < 2e-3)\ncheck(\"P(y = 1) = 1 - P(y = 0)\", np.isclose(np.mean(y == 1), 1 - f0))\n\npdf = lambda t: np.exp(-t ** 2 / 2) / np.sqrt(2 * np.pi)\nx = rng.standard_normal(10**7)\na = 0.5\nrel_errs = []\nfor width in (0.5, 0.1, 0.02):\n    p_emp = np.mean((x >= a) & (x <= a + width))\n    rel_errs.append(abs(p_emp - pdf(a) * width) / p_emp)\nprint(f\"    relative approximation error for |b-a|=0.5,0.1,0.02: \"\n      f\"{rel_errs[0]:.3f}, {rel_errs[1]:.3f}, {rel_errs[2]:.3f}\")\ncheck(\"the pdf approximation improves as the interval shrinks\",\n      rel_errs[0] > rel_errs[1] > rel_errs[2])\ncheck(\"relative error below 1% for |b - a| = 0.02\", rel_errs[2] < 0.01)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-picture]** The distribution visualized: the joint pdf of a data point (feature, label) as a grayscale over the plane, a two-component Gaussian mixture. Paralleling the entry's Fig. 1, three data sets of realizations of iid RVs with this distribution are drawn, and three hypothesis maps are learned from them by the same polynomial regression method. The checks verify that the plotted window carries the probability mass, that the data sets concentrate where the pdf is large, that the three learned maps differ, and that they differ most where the pdf is small."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-picture] grayscale pdf; three data sets, three learned maps\")\n# The joint pdf of a data point z = (x, y): the feature x follows a\n# two-component Gaussian mixture, and the label y is the value of a fixed\n# curve at x plus Gaussian noise.\nMIX_W = np.array([0.5, 0.5])\nMIX_MU = np.array([-2.0, 2.0])\nMIX_S = np.array([0.7, 0.9])\nNOISE = 0.35\nM_TRAIN, NR_SETS = 40, 3\n\n\ndef gauss_pdf(t, mu, s):\n    return np.exp(-((t - mu) ** 2) / (2 * s ** 2)) / (s * np.sqrt(2 * np.pi))\n\n\ndef pdf_feature(t):\n    return sum(w * gauss_pdf(t, mu, s) for w, mu, s in zip(MIX_W, MIX_MU, MIX_S))\n\n\ndef pdf_joint(t, u):\n    return pdf_feature(t) * gauss_pdf(u, np.tanh(t), NOISE)\n\n\ndef draw_points(m):\n    comp = rng.choice(2, size=m, p=MIX_W)\n    xx = rng.normal(MIX_MU[comp], MIX_S[comp])\n    yy = np.tanh(xx) + NOISE * rng.standard_normal(m)\n    return xx, yy\n\n\n# three data sets of iid draws, and the map the same polynomial regression\n# method learns from each (the degree-3 polynomial minimizing the average\n# squared deviation on the data set, available in closed form)\nx_curve = np.linspace(-4.4, 4.4, 200)\nsets, curves = [], []\nfor _ in range(NR_SETS):\n    x_tr, y_tr = draw_points(M_TRAIN)\n    design = np.vander(x_tr, 4)\n    w_hat = np.linalg.solve(design.T @ design, design.T @ y_tr)\n    sets.append((x_tr, y_tr))\n    curves.append(np.polyval(w_hat, x_curve))\ncurves = np.array(curves)\n\nxs = np.linspace(-4.4, 4.4, 61)\nys = np.linspace(-1.9, 1.9, 41)\ngrid_x, grid_y = np.meshgrid(xs, ys)\ngrid_p = pdf_joint(grid_x, grid_y)\n\nmass = float(grid_p.sum() * (xs[1] - xs[0]) * (ys[1] - ys[0]))\nprint(f\"    probability mass inside the plotted window: {mass:.4f}\")\ncheck(\"the plotted window carries the probability mass (within 2%)\",\n      abs(mass - 1.0) < 0.02)\n\nfor i, (x_tr, y_tr) in enumerate(sets, 1):\n    check(f\"data set {i} concentrates where the pdf is large\",\n          float(pdf_joint(x_tr, y_tr).mean()) > float(grid_p.mean()))\n\ngap01 = float(np.max(np.abs(curves[0] - curves[1])))\ncheck(\"the three learned maps differ\", gap01 > 0.05)\nspread_curve = curves.std(axis=0)                  # pointwise over the maps\ndens = pdf_feature(x_curve)\nlow, high = dens < np.median(dens), dens >= np.median(dens)\nprint(f\"    map spread: {spread_curve[low].mean():.3f} where the pdf is \"\n      f\"small, {spread_curve[high].mean():.3f} where it is large\")\ncheck(\"the maps differ most where the pdf is small\",\n      spread_curve[low].mean() > spread_curve[high].mean())\n\nnp.savetxt(OUT_DIR / \"probdist_density.csv\",\n           np.column_stack([grid_x.ravel(), grid_y.ravel(), grid_p.ravel()]),\n           delimiter=\",\", header=\"x,y,p\", comments=\"\", fmt=\"%.6f\")\nfor i, (x_tr, y_tr) in enumerate(sets, 1):\n    np.savetxt(OUT_DIR / f\"probdist_train{i}.csv\",\n               np.column_stack([x_tr, y_tr]),\n               delimiter=\",\", header=\"x,y\", comments=\"\", fmt=\"%.6f\")\nnp.savetxt(OUT_DIR / \"probdist_hypotheses.csv\",\n           np.column_stack([x_curve, curves[0], curves[1], curves[2]]),\n           delimiter=\",\", header=\"x,yhat1,yhat2,yhat3\", comments=\"\",\n           fmt=\"%.6f\")\n\n# ------------------------------------------------------------ preview\nfig, ax = plt.subplots(1, 3, figsize=(12.6, 3.0))\nfor m, style in ((25, \"--\"), (400, \"-\")):\n    ax[0].hist(outputs[m], bins=60, histtype=\"step\", density=True,\n               color=\"k\", linestyle=style, label=f\"m = {m}\")\nax[0].set_xlabel(\"output of the method (the average)\")\nax[0].set_ylabel(\"density over data sets\")\nax[0].set_title(\"outputs of one method over 4000 data sets\")\nax[0].legend(frameon=False)\n\nt = np.linspace(-4, 4, 400)\nax[1].plot(t, pdf(t), \"k-\")\nax[1].fill_between(t, pdf(t), where=(t >= a) & (t <= a + 0.5),\n                   facecolor=\"none\", hatch=\"///\", edgecolor=\"k\")\nax[1].set_xlabel(\"value of the RV\")\nax[1].set_ylabel(\"probability density p\")\nax[1].set_title(\"P(x in [a, b]) is the shaded area\")\n\nax[2].pcolormesh(xs, ys, grid_p, cmap=\"Greys\", shading=\"nearest\",\n                 vmin=0.0, vmax=float(grid_p.max()) * 1.3)\nfor i, (ls, mk) in enumerate(zip((\"-\", \"--\", \":\"), (\"o\", \"s\", \"^\"))):\n    ax[2].plot(x_curve, curves[i], \"k\", linestyle=ls,\n               label=f\"learned map {i + 1}\")\n    ax[2].scatter(sets[i][0], sets[i][1], s=18, marker=mk,\n                  facecolors=(\"black\", \"white\", \"0.5\")[i],\n                  edgecolors=(\"white\", \"black\", \"black\")[i],\n                  linewidths=0.4, label=f\"data set {i + 1}\", zorder=3)\nax[2].set_xlabel(\"feature x\")\nax[2].set_ylabel(\"label y\")\nax[2].set_title(\"pdf (grayscale), three data sets, three learned maps\")\nax[2].legend(frameon=False, loc=\"upper left\", fontsize=7)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"probdist.png\", dpi=110)\n\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}