{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "randomexperiment.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# random experiment \u2014 Python demo\n\nNumerical companion to the entry [random experiment](https://dictionaryofml.org/terms/randomexperiment.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 claims with three random experiments: rolling a die, whose probability measure is a model that the relative frequencies of many rolls check; the weather of a day at Krems an der Donau, whose outcome is read through a random variable (RV) whose probability distribution the probability measure on the sample space induces; and the repetition of an experiment, whose outcomes are sequences, whose measurements are i.i.d. RVs, and whose long-run behavior the law of large numbers (LLN) and the central limit theorem (CLT) make precise. 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/randomexperiment.py`](https://dictionaryofml.org/terms/randomexperiment.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(), \"randomexperiment.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nrandomexperiment.py \u2014 numerical companion to the glossary entry 'random\nexperiment'.\n\nPurpose\n-------\nBacks the entry's claims with three random experiments: rolling a die,\nwhose probability measure is a model that the relative frequencies of\nmany rolls check; the weather of a day at Krems an der Donau, whose\noutcome is read through a random variable (RV) whose probability\ndistribution the probability measure on the sample space induces; and\nthe repetition of an experiment, whose outcomes are sequences, whose\nmeasurements are i.i.d. RVs, and whose long-run behavior the law of\nlarge numbers (LLN) and the central limit theorem (CLT) make precise.\nSelf-contained (numpy/matplotlib only), fixed seed.\n\nSetup\n-----\nDie: sample space {1, ..., 6}; the uniform probability measure is the\nmodel of a fair die, a second measure with P(6) = 1/4 the model of a\nloaded one.  Weather: the sample space is the set of 365 days of 2023 at\nKrems with the uniform measure, the RV x the maximum temperature of a\nday (minimum and maximum downloaded from GeoSphere Austria); the days of\nJanuary 2023 are the entry's opening figure.  Repetition: m rolls of the\nfair die, m = 1 to 5000.\n\nBlocks\n------\n[B-die]     The probability measure is a property of the model, not of\n            the process: 6000 rolls of a fair die and of a loaded one\n            produce relative frequencies of a six near 1/6 and 1/4, and\n            the frequencies tell the two models apart.\n[B-january] The weather of a day at Krems as a random experiment, one\n            experiment per day: the minimum and the maximum temperature of\n            the 31 days of January 2023 are its measurements, the stems of\n            the entry's opening figure.  Read the other way, the whole\n            month is one random experiment whose outcome is the sequence.\n[B-rv]      The outcome of the weather experiment is a day; the\n            measurement x(omega) is its maximum temperature.  The induced\n            probability distribution P^(x)(B) = P(x in B) is the fraction\n            of days with x in B, checked for B = (-inf, 20].\n[B-repeat]  Repeating the die m times: the outcomes are sequences, the\n            product sample space has 6^m elements, the measurements are\n            i.i.d. RVs.  The average of the m rolls approaches 3.5 (LLN),\n            and the standardized sum of 30 rolls follows the standard\n            Gaussian bell (CLT), to within 0.01 in the cdf.\n[B-ml]      Two random experiments inside ML methods: the draw of a\n            subset of the training set for one SGD step, whose subset\n            average of the loss slopes has the full average as its\n            expectation, and the bootstrap resample of a\n            training set, in which a given data point appears with\n            probability near 1 - 1/e.\n\nOutputs\n-------\nrandomexperiment_january.csv : day of January 2023 at Krems, minimum and\n                            maximum temperature\nrandomexperiment_die.csv  : number of rolls, running relative frequency of\n                            a six for the fair and the loaded die\nrandomexperiment_lln.csv  : m, average of m rolls\nrandomexperiment_clt.csv  : standardized sum of 30 rolls, its cdf over\n                            10000 repetitions, and the Gaussian cdf\nrandomexperiment.png      : matplotlib preview (checking only)\n\"\"\"\n\nimport json\nimport math\nimport urllib.request\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(2020)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-die]** The probability measure is a property of the model, not of the process: 6000 rolls of a fair die and of a loaded one produce relative frequencies of a six near 1/6 and 1/4, and the frequencies tell the two models apart."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "N_ROLLS = 6000\nOMEGA = np.arange(1, 7)                              # the sample space of one roll\nP_FAIR = np.full(6, 1 / 6)                           # the model of a fair die\nP_LOADED = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.25])\nrolls_fair = rng.choice(OMEGA, N_ROLLS, p=P_FAIR)\nrolls_loaded = rng.choice(OMEGA, N_ROLLS, p=P_LOADED)\nfreq_fair = np.cumsum(rolls_fair == 6) / np.arange(1, N_ROLLS + 1)\nfreq_loaded = np.cumsum(rolls_loaded == 6) / np.arange(1, N_ROLLS + 1)\ncheck(f\"[B-die]     after {N_ROLLS} rolls the relative frequency of a six is {freq_fair[-1]:.3f} \"\n      f\"for the fair die (model 1/6) and {freq_loaded[-1]:.3f} for the loaded one (model 1/4)\",\n      abs(freq_fair[-1] - 1 / 6) < 0.02 and abs(freq_loaded[-1] - 0.25) < 0.02)\ncheck(\"[B-die]     the frequencies tell the two models apart: the fair die's frequency is \"\n      \"closer to 1/6 than to 1/4, the loaded die's the other way round\",\n      abs(freq_fair[-1] - 1 / 6) < abs(freq_fair[-1] - 0.25)\n      and abs(freq_loaded[-1] - 0.25) < abs(freq_loaded[-1] - 1 / 6))\nwith open(OUT_DIR / \"randomexperiment_die.csv\", \"w\") as fh:\n    fh.write(\"rolls,fair,loaded\\n\")\n    for n in list(range(1, 100)) + list(range(100, N_ROLLS + 1, 20)):\n        fh.write(f\"{n},{freq_fair[n - 1]:.4f},{freq_loaded[n - 1]:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-january]** The weather of a day at Krems as a random experiment, one experiment per day: the minimum and the maximum temperature of the 31 days of January 2023 are its measurements, the stems of the entry's opening figure. Read the other way, the whole month is one random experiment whose outcome is the sequence."
  },
  {
   "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&start=2023-01-01&end=2023-12-31\")\nwith urllib.request.urlopen(url, timeout=120) as resp:\n    rec = json.load(resp)\npar = rec[\"features\"][0][\"properties\"][\"parameters\"]\ntmin = np.array(par[\"tlmin\"][\"data\"], dtype=float)\ntmax = np.array(par[\"tlmax\"][\"data\"], dtype=float)\ndays = np.array(rec[\"timestamps\"])\njan = np.array([d.startswith(\"2023-01\") for d in days])\ncheck(f\"[B-january] {int(jan.sum())} days of January 2023 at Krems, each the outcome of one random \"\n      \"experiment, measured by its minimum and its maximum temperature\",\n      jan.sum() == 31 and not np.isnan(tmin[jan]).any() and not np.isnan(tmax[jan]).any())\ncheck(f\"[B-january] the two measurements of a day are ordered, minimum below maximum, on every day \"\n      f\"(smallest gap {np.min(tmax[jan] - tmin[jan]):.1f} degrees)\", bool(np.all(tmax[jan] > tmin[jan])))\nwith open(OUT_DIR / \"randomexperiment_january.csv\", \"w\") as fh:\n    fh.write(\"day,tmin,tmax\\n\")\n    for k, (a, b) in enumerate(zip(tmin[jan], tmax[jan]), 1):\n        fh.write(f\"{k},{a:.1f},{b:.1f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-rv]** The outcome of the weather experiment is a day; the measurement x(omega) is its maximum temperature. The induced probability distribution P^(x)(B) = P(x in B) is the fraction of days with x in B, checked for B = (-inf, 20]."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "check(f\"[B-rv]      the sample space of the weather experiment: {len(days)} days of 2023 at Krems, \"\n      \"each an outcome omega with the uniform probability measure\",\n      len(days) == 365 and not np.isnan(tmax).any())\nB = 20.0\ninduced = float(np.mean(tmax <= B))                  # P^(x)((-inf, B]) = P({omega : x(omega) <= B})\ncount = int((tmax <= B).sum())\ncheck(f\"[B-rv]      the RV x(omega) = maximum temperature of the day induces the distribution \"\n      f\"P^(x)((-inf, {B:g}]) = P(x <= {B:g}) = {count}/365 = {induced:.3f}\", abs(induced - count / 365) < 1e-12)\ncheck(f\"[B-rv]      the realization on 2023-07-01 is x = {tmax[181]:.1f} degrees, one value of the RV\",\n      days[181].startswith(\"2023-07-01\"))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-repeat]** Repeating the die m times: the outcomes are sequences, the product sample space has 6^m elements, the measurements are i.i.d. RVs. The average of the m rolls approaches 3.5 (LLN), and the standardized sum of 30 rolls follows the standard Gaussian bell (CLT), to within 0.01 in the cdf."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "M_MAX = 5000\nseq = rng.choice(OMEGA, M_MAX, p=P_FAIR)             # one outcome of the repeated experiment\navg = np.cumsum(seq) / np.arange(1, M_MAX + 1)\ncheck(f\"[B-repeat]  repeating the roll m times has 6^m outcomes (6^30 = {6 ** 30:.1e}); the average \"\n      f\"of m rolls approaches 3.5: {avg[9]:.2f} at m = 10, {avg[99]:.2f} at m = 100, \"\n      f\"{avg[-1]:.3f} at m = {M_MAX} (LLN)\", abs(avg[-1] - 3.5) < 0.05)\nM_CLT, REPS = 30, 10000\nsums = rng.choice(OMEGA, (REPS, M_CLT), p=P_FAIR).sum(1)\nmu, sd = 3.5, math.sqrt(M_CLT * 35 / 12)\n# the sum takes integer values; its cdf at the integer s is compared with the\n# Gaussian cdf at the standardized value of s + 1/2, the upper edge of that lattice cell\nsvals = np.arange(int(M_CLT * mu - 3 * sd), int(M_CLT * mu + 3 * sd) + 1)\ngrid = (svals + 0.5 - M_CLT * mu) / sd\nemp_cdf = np.array([np.mean(sums <= sv) for sv in svals])\ngauss_cdf = np.array([0.5 * (1 + math.erf(g / math.sqrt(2))) for g in grid])\ngap = float(np.max(np.abs(emp_cdf - gauss_cdf)))\ncheck(f\"[B-repeat]  the standardized sum of {M_CLT} rolls follows the standard Gaussian bell: largest \"\n      f\"gap between the cdfs over {REPS} repetitions {gap:.3f} (CLT)\", gap < 0.01)\nwith open(OUT_DIR / \"randomexperiment_lln.csv\", \"w\") as fh:\n    fh.write(\"m,average\\n\")\n    for m in list(range(1, 100)) + list(range(100, M_MAX + 1, 20)):\n        fh.write(f\"{m},{avg[m - 1]:.4f}\\n\")\nwith open(OUT_DIR / \"randomexperiment_clt.csv\", \"w\") as fh:\n    fh.write(\"z,cdf_sum,cdf_gauss\\n\")\n    for g, a, b in zip(grid, emp_cdf, gauss_cdf):\n        fh.write(f\"{g:.3f},{a:.4f},{b:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-ml]** Two random experiments inside ML methods: the draw of a subset of the training set for one SGD step, whose subset average of the loss slopes has the full average as its expectation, and the bootstrap resample of a training set, in which a given data point appears with probability near 1 - 1/e."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "m_train, n_sub = 365, 20\nx_tr = tmax[:-1]; y_tr = tmax[1:]                    # predict tomorrow's maximum from today's\nw = np.array([0.0, 0.5])\nfull_grad = 2 * np.mean((w[0] + w[1] * x_tr - y_tr)[:, None] * np.c_[np.ones(364), x_tr], 0)\nsub_slopes = np.array([2 * np.mean((w[0] + w[1] * x_tr[i] - y_tr[i])[:, None] * np.c_[np.ones(n_sub), x_tr[i]], 0)\n                       for i in (rng.choice(364, n_sub, replace=False) for _ in range(4000))])\ncheck(f\"[B-ml]      the subset draw of SGD is a random experiment with C(364, {n_sub}) outcomes; the \"\n      f\"subset average of the loss slopes over 4000 draws, {np.round(sub_slopes.mean(0), 2)}, has the \"\n      f\"full average {np.round(full_grad, 2)} as its expectation\",\n      np.allclose(sub_slopes.mean(0), full_grad, rtol=0.05, atol=0.05))\ninc = np.mean([len(set(rng.integers(0, m_train, m_train))) / m_train for _ in range(2000)])\ncheck(f\"[B-ml]      a bootstrap resample of the {m_train} days holds a given day with probability \"\n      f\"{inc:.3f}, near 1 - 1/e = {1 - 1 / math.e:.3f}\", abs(inc - (1 - 1 / math.e)) < 0.01)\n\n# -------------------------------------------------------------- preview\nfig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12.5, 3.7))\nn = np.arange(1, N_ROLLS + 1)\nax1.plot(n, freq_fair, \"k-\", label=\"fair die\"); ax1.plot(n, freq_loaded, \"k--\", label=\"loaded die\")\nax1.axhline(1 / 6, color=\"0.6\", lw=0.8, ls=\":\"); ax1.axhline(0.25, color=\"0.6\", lw=0.8, ls=\":\")\nax1.set_xscale(\"log\"); ax1.set_xlabel(\"number of rolls\"); ax1.set_ylabel(\"relative frequency of a six\")\nax1.set_title(\"relative frequencies check the model\"); ax1.legend(frameon=False, fontsize=8)\ndj = np.arange(1, 32)\nax2.vlines(dj, 0, tmax[jan], colors=\"k\", lw=1.0)\nax2.vlines(dj, 0, tmin[jan], colors=\"k\", lw=1.0, linestyles=\"dotted\")\nax2.plot(dj, tmax[jan], \"ko\", ms=4, label=\"maximum temperature (solid stems)\")\nax2.plot(dj, tmin[jan], \"ks\", ms=4, mfc=\"white\", label=\"minimum temperature (dotted stems)\")\nax2.axhline(0, color=\"0.6\", lw=0.8)\nax2.set_xlabel(\"day of January 2023\"); ax2.set_ylabel(\"temperature at Krems (deg C)\")\nax2.set_title(\"one random experiment per day, two measurements\"); ax2.legend(frameon=False, fontsize=8)\nax3.plot(grid, emp_cdf, \"k-\", label=f\"standardized sum of {M_CLT} rolls\")\nax3.plot(grid, gauss_cdf, \"k--\", label=\"standard Gaussian\")\nax3.set_xlabel(\"$z$\"); ax3.set_ylabel(\"cdf\"); ax3.set_title(\"repetition: the CLT\"); ax3.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"randomexperiment.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 / 'randomexperiment_january.csv'}, {OUT_DIR / 'randomexperiment_die.csv'}, {OUT_DIR / 'randomexperiment_lln.csv'}, \"\n      f\"{OUT_DIR / 'randomexperiment_clt.csv'}, {OUT_DIR / 'randomexperiment.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}