{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "em.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# expectation\u2013maximization \u2014 Python demo\n\nNumerical companion to the entry [expectation\u2013maximization](https://dictionaryofml.org/terms/em.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 'em' (expectation-maximization). The GeoSphere Austria weather station Krems (station id 3805) records the minimum air temperature of each night; this script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them to em_temps.csv. The nightly minima scatter around two regimes -- cold-season and warm-season nights -- and a Gaussian mixture model (GMM) with two components captures exactly such a distribution. Maximizing its likelihood has no closed-form solution, so the model parameters are fitted by the EM algorithm.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/em.py`](https://dictionaryofml.org/terms/em.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(), \"em.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"EM fits a two-component GMM to a year of nightly temperatures at\nKrems: the negative log-likelihood never increases, iteration by\niteration.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'em'\n(expectation-maximization).  The GeoSphere Austria weather station\nKrems (station id 3805) records the minimum air temperature of each\nnight; this script downloads the records for 2024 from the GeoSphere\ndata hub (dataset klima-v2-1d) and writes them to em_temps.csv.  The\nnightly minima scatter around two regimes -- cold-season and\nwarm-season nights -- and a Gaussian mixture model (GMM) with two\ncomponents captures exactly such a distribution.  Maximizing its\nlikelihood has no closed-form solution, so the model parameters are\nfitted by the EM algorithm.\n\nThe demo checks the entry's central claims: each EM iteration\nminimizes a surrogate objective that upper-bounds the negative\nlog-likelihood and is tight at the current iterate, so the negative\nlog-likelihood never increases; the iteration stops at a fixed point,\nwhere the parameters minimize their own surrogate.\n\nDeterministic: the initialization is the 25th/75th percentile of the\ndata (no randomness).  Self-contained: numpy + matplotlib only\n(stdlib urllib for the download).\n\nBlocks\n------\n[B-fetch] Download the 366 nightly minimum temperatures at Krems for\n          2024 and write them to em_temps.csv; check the count and one\n          pinned value against the archive.\n[B-em]    Run EM for a two-component GMM: E-step (posterior\n          probabilities of the two components), M-step (re-weighted\n          means, variances, and component probabilities).  Check that\n          the negative log-likelihood never increases and that the\n          iteration reaches a fixed point.\n[B-fit]   The fitted mixture: two well-separated component means (a\n          cold-season and a warm-season regime); write the histogram,\n          the fitted densities, and the negative log-likelihood per\n          iteration for the entry's figure.\n\nOutputs\n-------\nem_temps.csv   : date, tmin for the 366 nights of 2024\nem_hist.csv    : t, freq -- normalized histogram of the temperatures\nem_density.csv : t, mix, comp1, comp2 -- fitted mixture and components\nem_loglik.csv  : iter, nll -- negative log-likelihood per EM iteration\nem.png         : preview (checking only) -- histogram with the fitted\n                 densities, and the monotone negative log-likelihood\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 nightly minimum temperatures at Krems for 2024 and write them to em_temps.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&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)\nstamps = [t[:10] for t in payload[\"timestamps\"]]\ntmin = np.array(payload[\"features\"][0][\"properties\"][\"parameters\"]\n                [\"tlmin\"][\"data\"], dtype=float)\nwith open(OUT_DIR / \"em_temps.csv\", \"w\") as f:\n    f.write(\"date,tmin\\n\")\n    for day, t in zip(stamps, tmin):\n        f.write(f\"{day},{t}\\n\")\ncheck(\"[B-fetch] 366 nightly minima downloaded for 2024\",\n      len(tmin) == 366)\ncheck(\"[B-fetch] the record matches the archive (Feb 1: -3.8)\",\n      stamps[31] == \"2024-02-01\" and abs(tmin[31] - (-3.8)) < 1e-9)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-em]** Run EM for a two-component GMM: E-step (posterior probabilities of the two components), M-step (re-weighted means, variances, and component probabilities). Check that the negative log-likelihood never increases and that the iteration reaches a fixed point."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def normal_pdf(t, mean, var):\n    return np.exp(-0.5 * (t - mean) ** 2 / var) / np.sqrt(2 * np.pi * var)\n\n\ndef neg_log_likelihood(t, p, means, variances):\n    mix = sum(p[c] * normal_pdf(t, means[c], variances[c]) for c in (0, 1))\n    return float(-np.log(mix).sum())\n\n\nmeans = np.percentile(tmin, [25.0, 75.0])          # deterministic start\nvariances = np.array([tmin.var(), tmin.var()])\np = np.array([0.5, 0.5])\n\nnll_trace = [neg_log_likelihood(tmin, p, means, variances)]\nfor _ in range(200):\n    # E-step: posterior probability of each component per night\n    joint = np.stack([p[c] * normal_pdf(tmin, means[c], variances[c])\n                      for c in (0, 1)])\n    posterior = joint / joint.sum(axis=0)\n    # M-step: re-weighted component probabilities, means, variances\n    weight = posterior.sum(axis=1)\n    p = weight / len(tmin)\n    means = (posterior * tmin).sum(axis=1) / weight\n    variances = (posterior * (tmin - means[:, None]) ** 2).sum(axis=1) / weight\n    nll_trace.append(neg_log_likelihood(tmin, p, means, variances))\nnll_trace = np.array(nll_trace)\n\ncheck(\"[B-em] the negative log-likelihood never increases\",\n      bool(np.all(np.diff(nll_trace) <= 1e-9)))\ncheck(\"[B-em] the iteration reaches a fixed point (last update tiny)\",\n      abs(nll_trace[-1] - nll_trace[-2]) < 1e-10)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fit]** The fitted mixture: two well-separated component means (a cold-season and a warm-season regime); write the histogram, the fitted densities, and the negative log-likelihood per iteration for the entry's figure."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "order = np.argsort(means)\np, means, variances = p[order], means[order], variances[order]\nprint(f\"  cold regime: mean {means[0]:+.1f} C (probability {p[0]:.2f}), \"\n      f\"warm regime: mean {means[1]:+.1f} C (probability {p[1]:.2f})\")\ncheck(\"[B-fit] the two component means are well separated\",\n      means[1] - means[0] > 5.0)\n\ncounts, edges = np.histogram(tmin, bins=24, density=True)\ncenters = 0.5 * (edges[:-1] + edges[1:])\nnp.savetxt(OUT_DIR / \"em_hist.csv\",\n           np.stack([centers, counts], 1), delimiter=\",\",\n           header=\"t,freq\", comments=\"\", fmt=\"%.4f\")\ngrid = np.linspace(tmin.min() - 2, tmin.max() + 2, 300)\ncomp = [p[c] * normal_pdf(grid, means[c], variances[c]) for c in (0, 1)]\nnp.savetxt(OUT_DIR / \"em_density.csv\",\n           np.stack([grid, comp[0] + comp[1], comp[0], comp[1]], 1),\n           delimiter=\",\", header=\"t,mix,comp1,comp2\", comments=\"\",\n           fmt=\"%.5f\")\nnp.savetxt(OUT_DIR / \"em_loglik.csv\",\n           np.stack([np.arange(len(nll_trace)), nll_trace], 1),\n           delimiter=\",\", header=\"iter,nll\", comments=\"\", fmt=\"%.4f\")\n\nfig, axes = plt.subplots(1, 2, figsize=(9, 3.4))\nax = axes[0]\nax.bar(centers, counts, width=edges[1] - edges[0], color=\"0.8\",\n       edgecolor=\"0.5\", label=\"nightly minima 2024\")\nax.plot(grid, comp[0] + comp[1], \"k-\", label=\"fitted GMM\")\nax.plot(grid, comp[0], \"k--\", label=\"cold-season component\")\nax.plot(grid, comp[1], \"k:\", label=\"warm-season component\")\nax.set_xlabel(\"nightly minimum temperature (\u00b0C)\")\nax.set_ylabel(\"relative frequency\")\nax.set_title(\"two-component GMM fitted by EM (Krems, 2024)\")\nax.legend(frameon=False, fontsize=8)\nax = axes[1]\nax.plot(np.arange(len(nll_trace)), nll_trace, \"k-\")\nax.set_xlabel(\"EM iteration\")\nax.set_ylabel(\"negative log-likelihood\")\nax.set_title(\"monotone descent to a fixed point\")\nax.set_xlim(0, 30)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"em.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 \"\"))"
  }
 ]
}