B-majorize writes when the script runsDictionary of Applied Machine Learning · expectation–maximization (EM)
Numerical companion to the entry expectation–maximization (EM): 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 without installing anything:uv run https://dictionaryofml.org/terms/em.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download em.py into a folder and run uv run em.py there. With NumPy and Matplotlib already installed, python3 em.py, from any directory — it writes its output files into the current directory. 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.
[B-majorize] The majorize-minimize view on one parameter: with the
cluster probabilities and the warm-season component fixed at
their learned values and the cold-season variance fixed at a
smaller value (so that the surrogates differ visibly from the
risk), the empirical risk as a function of the cold-season
mean, and the surrogate built by the E-step at each of the
first iterates of EM on that mean. Check that each surrogate
lies above the empirical risk, touches it at the current
iterate, and is minimized at the next iterate.
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_majorant_risk.csv : w, risk -- empirical risk (average negative
log-likelihood) as a function of the cold-season mean w
em_majorant_surrogate1.csv, _surrogate2.csv, _surrogate3.csv : w, q --
the surrogate built at the first, second and third
iterate of EM on the cold-season mean
em_majorant_iterates.csv : iter, w, risk -- the iterates of EM on the
cold-season mean and their empirical risk
em.png : preview (checking only) -- histogram with the fitted
densities, the monotone negative log-likelihood, and
the surrogates of the majorize-minimize view
"""
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, variance 39.7 (probability 0.83); warm regime: mean +17.2 C, variance 3.3 (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, variance {variances[0]:.1f} "
f"(probability {p[0]:.2f}); warm regime: mean {means[1]:+.1f} C, "
f"variance {variances[1]:.1f} (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")
[ok] [B-fit] the two component means are well separated EM on the cold-season mean (variance fixed at 16) from w = -8.00, +2.27, +3.96, +4.18; minimizer of the risk on the grid +4.20
The majorize-minimize view on one parameter: with the cluster probabilities and the warm-season component fixed at their learned values and the cold-season variance fixed at a smaller value (so that the surrogates differ visibly from the risk), the empirical risk as a function of the cold-season mean, and the surrogate built by the E-step at each of the first iterates of EM on that mean. Check that each surrogate lies above the empirical risk, touches it at the current iterate, and is minimized at the next iterate.
VAR_COLD = 16.0 # cold-season variance for the illustration (learned: ~40)
def risk_of_cold_mean(w):
"""Average negative log-likelihood with the cold-season mean at w."""
mix = (p[0] * normal_pdf(tmin, w, VAR_COLD)
+ p[1] * normal_pdf(tmin, means[1], variances[1]))
return float(-np.log(mix).mean())
def surrogate_of_cold_mean(w, w_cur):
"""E-step at w_cur, then the surrogate Q(w | w_cur) on the grid w."""
cold = p[0] * normal_pdf(tmin, w_cur, VAR_COLD)
warm = p[1] * normal_pdf(tmin, means[1], variances[1])
u = cold / (cold + warm) # posterior of the cold regime
entropy = xlogx(u) + xlogx(1.0 - u) # sum_y u_y log u_y per night
q = np.empty_like(w)
for i, wi in enumerate(w):
cold_w = p[0] * normal_pdf(tmin, wi, VAR_COLD)
q[i] = -np.mean(u * np.log(cold_w) + (1.0 - u) * np.log(warm)
- entropy)
return q, float((u * tmin).sum() / u.sum())
def xlogx(a):
"""a log a with the convention 0 log 0 = 0."""
out = np.zeros_like(a)
pos = a > 0
out[pos] = a[pos] * np.log(a[pos])
return out
W_START = -8.0 # a poor first guess
w_grid = np.linspace(-13.0, 15.0, 561)
risk_grid = np.array([risk_of_cold_mean(w) for w in w_grid])
np.savetxt(OUT_DIR / "em_majorant_risk.csv",
np.stack([w_grid, risk_grid], 1), delimiter=",",
header="w,risk", comments="", fmt="%.5f")
iterates = [W_START]
surrogates = []
for i in (1, 2, 3):
w_cur = iterates[-1]
q, w_next = surrogate_of_cold_mean(w_grid, w_cur)
surrogates.append(q)
iterates.append(w_next)
np.savetxt(OUT_DIR / f"em_majorant_surrogate{i}.csv",
np.stack([w_grid, q], 1), delimiter=",",
header="w,q", comments="", fmt="%.5f")
iter_risk = np.array([risk_of_cold_mean(w) for w in iterates])
np.savetxt(OUT_DIR / "em_majorant_iterates.csv",
np.stack([np.arange(len(iterates)), iterates, iter_risk], 1),
delimiter=",", header="iter,w,risk", comments="", fmt="%.5f")
print(" EM on the cold-season mean (variance fixed at "
f"{VAR_COLD:g}) from w = "
+ ", ".join(f"{w:+.2f}" for w in iterates)
+ f"; minimizer of the risk on the grid {w_grid[np.argmin(risk_grid)]:+.2f}")
check("[B-majorize] each surrogate lies above the empirical risk",
all(bool(np.all(q >= risk_grid - 1e-9)) for q in surrogates))
check("[B-majorize] each surrogate touches the risk at its iterate",
all(abs(np.interp(iterates[i], w_grid, surrogates[i])
- risk_of_cold_mean(iterates[i])) < 2e-4 for i in range(3)))
check("[B-majorize] each surrogate is minimized at the next iterate",
all(abs(w_grid[np.argmin(surrogates[i])] - iterates[i + 1]) < 0.06
for i in range(3)))
check("[B-majorize] the empirical risk decreases along the iterates",
bool(np.all(np.diff(iter_risk) < 0)))
fig, axes = plt.subplots(1, 3, figsize=(13, 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, "ko", markersize=3,
linestyle="none")
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)
ax = axes[2]
ax.plot(w_grid, risk_grid, "k-", label="empirical risk")
for i, q in enumerate(surrogates):
ax.plot(w_grid, q, "k--", linewidth=0.9,
label="surrogate at an iterate" if i == 0 else None)
ax.plot(iterates, iter_risk, "ko", markersize=4, linestyle="none",
label="iterates")
ax.set_ylim(risk_grid.min() - 0.1, risk_grid.max() + 1.0)
ax.set_xlabel("cold-season mean w (°C)")
ax.set_ylabel("empirical risk")
ax.set_title("surrogates majorize the risk")
ax.legend(frameon=False, fontsize=8)
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-majorize] each surrogate lies above the empirical risk [ok] [B-majorize] each surrogate touches the risk at its iterate [ok] [B-majorize] each surrogate is minimized at the next iterate [ok] [B-majorize] the empirical risk decreases along the iterates 9/9 checks passed
B-majorize writes when the script runs