B-fit writes when the script runsDictionary of Applied Machine Learning · expectation–maximization
Numerical companion to the entry expectation–maximization: it recomputes what the entry states and prints one line per check
Numerical 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.
Run it with python3 em.py, from any directory — it writes its output files into the current directory. Requires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Download em.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""EM fits a two-component GMM to a year of nightly temperatures at
Krems: the negative log-likelihood never increases, iteration by
iteration.
Purpose
-------
Numerical 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.
The demo checks the entry's central claims: each EM iteration
minimizes a surrogate objective that upper-bounds the negative
log-likelihood and is tight at the current iterate, so the negative
log-likelihood never increases; the iteration stops at a fixed point,
where the parameters minimize their own surrogate.
Deterministic: the initialization is the 25th/75th percentile of the
data (no randomness). Self-contained: numpy + matplotlib only
(stdlib urllib for the download).
Blocks
------
[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.
[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.
[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.
Outputs
-------
em_temps.csv : date, tmin for the 366 nights of 2024
em_hist.csv : t, freq -- normalized histogram of the temperatures
em_density.csv : t, mix, comp1, comp2 -- fitted mixture and components
em_loglik.csv : iter, nll -- negative log-likelihood per EM iteration
em.png : preview (checking only) -- histogram with the fitted
densities, and the monotone negative log-likelihood
"""
import json
import urllib.request
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
OUT_DIR = Path(__file__).parent
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
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.
URL = ("https://dataset.api.hub.geosphere.at/v1/station/historical/"
"klima-v2-1d?parameters=tlmin&station_ids=3805"
"&start=2024-01-01&end=2024-12-31")
with urllib.request.urlopen(URL, timeout=120) as resp:
payload = json.load(resp)
stamps = [t[:10] for t in payload["timestamps"]]
tmin = np.array(payload["features"][0]["properties"]["parameters"]
["tlmin"]["data"], dtype=float)
with open(OUT_DIR / "em_temps.csv", "w") as f:
f.write("date,tmin\n")
for day, t in zip(stamps, tmin):
f.write(f"{day},{t}\n")
check("[B-fetch] 366 nightly minima downloaded for 2024",
len(tmin) == 366)
check("[B-fetch] the record matches the archive (Feb 1: -3.8)",
stamps[31] == "2024-02-01" and abs(tmin[31] - (-3.8)) < 1e-9)
[ok] [B-fetch] 366 nightly minima downloaded for 2024 [ok] [B-fetch] the record matches the archive (Feb 1: -3.8)
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.
def normal_pdf(t, mean, var):
return np.exp(-0.5 * (t - mean) ** 2 / var) / np.sqrt(2 * np.pi * var)
def neg_log_likelihood(t, p, means, variances):
mix = sum(p[c] * normal_pdf(t, means[c], variances[c]) for c in (0, 1))
return float(-np.log(mix).sum())
means = np.percentile(tmin, [25.0, 75.0]) # deterministic start
variances = np.array([tmin.var(), tmin.var()])
p = np.array([0.5, 0.5])
nll_trace = [neg_log_likelihood(tmin, p, means, variances)]
for _ in range(200):
# E-step: posterior probability of each component per night
joint = np.stack([p[c] * normal_pdf(tmin, means[c], variances[c])
for c in (0, 1)])
posterior = joint / joint.sum(axis=0)
# M-step: re-weighted component probabilities, means, variances
weight = posterior.sum(axis=1)
p = weight / len(tmin)
means = (posterior * tmin).sum(axis=1) / weight
variances = (posterior * (tmin - means[:, None]) ** 2).sum(axis=1) / weight
nll_trace.append(neg_log_likelihood(tmin, p, means, variances))
nll_trace = np.array(nll_trace)
check("[B-em] the negative log-likelihood never increases",
bool(np.all(np.diff(nll_trace) <= 1e-9)))
check("[B-em] the iteration reaches a fixed point (last update tiny)",
abs(nll_trace[-1] - nll_trace[-2]) < 1e-10)
[ok] [B-em] the negative log-likelihood never increases [ok] [B-em] the iteration reaches a fixed point (last update tiny) cold regime: mean +5.1 C (probability 0.83), warm regime: mean +17.2 C (probability 0.17)
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.
order = np.argsort(means)
p, means, variances = p[order], means[order], variances[order]
print(f" cold regime: mean {means[0]:+.1f} C (probability {p[0]:.2f}), "
f"warm regime: mean {means[1]:+.1f} C (probability {p[1]:.2f})")
check("[B-fit] the two component means are well separated",
means[1] - means[0] > 5.0)
counts, edges = np.histogram(tmin, bins=24, density=True)
centers = 0.5 * (edges[:-1] + edges[1:])
np.savetxt(OUT_DIR / "em_hist.csv",
np.stack([centers, counts], 1), delimiter=",",
header="t,freq", comments="", fmt="%.4f")
grid = np.linspace(tmin.min() - 2, tmin.max() + 2, 300)
comp = [p[c] * normal_pdf(grid, means[c], variances[c]) for c in (0, 1)]
np.savetxt(OUT_DIR / "em_density.csv",
np.stack([grid, comp[0] + comp[1], comp[0], comp[1]], 1),
delimiter=",", header="t,mix,comp1,comp2", comments="",
fmt="%.5f")
np.savetxt(OUT_DIR / "em_loglik.csv",
np.stack([np.arange(len(nll_trace)), nll_trace], 1),
delimiter=",", header="iter,nll", comments="", fmt="%.4f")
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
ax = axes[0]
ax.bar(centers, counts, width=edges[1] - edges[0], color="0.8",
edgecolor="0.5", label="nightly minima 2024")
ax.plot(grid, comp[0] + comp[1], "k-", label="fitted GMM")
ax.plot(grid, comp[0], "k--", label="cold-season component")
ax.plot(grid, comp[1], "k:", label="warm-season component")
ax.set_xlabel("nightly minimum temperature (°C)")
ax.set_ylabel("relative frequency")
ax.set_title("two-component GMM fitted by EM (Krems, 2024)")
ax.legend(frameon=False, fontsize=8)
ax = axes[1]
ax.plot(np.arange(len(nll_trace)), nll_trace, "k-")
ax.set_xlabel("EM iteration")
ax.set_ylabel("negative log-likelihood")
ax.set_title("monotone descent to a fixed point")
ax.set_xlim(0, 30)
fig.tight_layout()
fig.savefig(OUT_DIR / "em.png", dpi=150)
failed = [name for name, ok in report if not ok]
print(f"{len(report) - len(failed)}/{len(report)} checks passed"
+ (f", FAILED: {failed}" if failed else ""))
[ok] [B-fit] the two component means are well separated 5/5 checks passed
B-fit writes when the script runs