B-plot writes when the script runsDictionary of Applied Machine Learning · Gaussian mixture model (GMM)
Updated on 2026-11-03
Numerical companion to the entry Gaussian mixture model (GMM): it recomputes what the entry states and prints one line per check
Numerical companion to the glossary entry 'gmm' (Gaussian mixture model). A photograph of the Oetscher massif (assets/oetscher.jpg, the same photograph the 'kmeans' entry uses) is cut into square patches. Each patch is a data point whose feature vector holds two numbers, how green and how blue the patch is on average, so the covariance matrix of a component is a genuine matrix rather than a single number and each component shows as an ellipse in the plane.
Run it without installing anything:uv run https://dictionaryofml.org/terms/gmm.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 gmm.py into a folder and run uv run gmm.py there. With NumPy and Matplotlib already installed, python3 gmm.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download gmm.py · Notebook · Open in Colab
Before running a downloaded copy, scan it: uvx bandit -q -t B102,B307,B602,B605 gmm.py (B102,B307,B602,B605 are the exec, eval and shell checks of bandit) reports nothing: the script runs no shell command and evaluates no generated code.
One cell per block of the script: the code, and what that code printed when it last ran here
"""A three-component GMM on the greenness and blueness of image patches:
meadow, foliage and sky, with graded membership where they meet.
Purpose
-------
Numerical companion to the glossary entry 'gmm' (Gaussian mixture
model). A photograph of the Oetscher massif (assets/oetscher.jpg, the
same photograph the 'kmeans' entry uses) is cut into square patches.
Each patch is a data point whose feature vector holds two numbers, how
green and how blue the patch is on average, so the covariance matrix of
a component is a genuine matrix rather than a single number and each
component shows as an ellipse in the plane.
The demo checks the entry's claims: the cluster probabilities sum to
one; the three components are the sunlit meadow, the dark foliage and
the sky with the mountain; each covariance matrix is symmetric positive semi-definite;
the posterior distribution grades the membership of a patch in each
cluster instead of assigning it to exactly one; and the EM algorithm
reduces to k-means when the cluster probabilities are equal and every
covariance matrix is a shrinking multiple of the identity matrix.
Deterministic: the components are initialized by splitting the patches
into three equal parts ordered by blueness (no randomness).
Self-contained: numpy + matplotlib only.
Blocks
------
[B-patches] Cut the photograph into 64x64 patches and measure how
green and how blue each one is; check the patch count.
[B-em] Fit the three-component GMM by the EM algorithm: E-step
(posterior probability of each cluster index per patch),
M-step (re-weighted cluster probabilities, means and
covariance matrices). Check that the cluster
probabilities sum to one and that the fit stops changing.
[B-components] The three components are the sunlit meadow, the dark
foliage, and the sky with the mountain. Check each against the average color
of its patches and against where they sit in the
photograph, and check that each covariance matrix is
symmetric positive semi-definite.
[B-soft] Soft clustering: the posterior distribution grades
membership. Check that most patches are graded
decisively and that the rest are shared between clusters.
[B-images] The photograph with the patch grid drawn on it, the
photograph at patch resolution, and one copy per
cluster whose brightness is that cluster's posterior
probability. Because the posterior probabilities of a
patch sum to one, the three copies add back up to the
photograph.
[B-kmeans] With equal cluster probabilities and covariance matrices
that shrink to a multiple of the identity matrix, the
posterior concentrates on the nearest mean and the EM
update becomes the k-means update, so the algorithm
reduces to k-means. Check both.
[B-sigma] The same constrained model for a small, a moderate and a
large common variance: the photograph dimmed by the
posterior probability of the vegetation cluster, next to
the k-means assignment to that cluster. Check that the
small variance reproduces the k-means assignment, the
moderate one grades the patches, and the large one
grades every patch equally.
[B-plot] Write the patch scatter and the three component ellipses
for the entry's figure, plus the preview.
Outputs
-------
gmm_patches.csv : green, blue -- every patch
gmm_points.csv : x1, x2 -- every patch
gmm_ellipse1.csv : x1, x2 -- contour of the meadow component
gmm_ellipse2.csv : x1, x2 -- contour of the foliage component
gmm_ellipse3.csv : x1, x2 -- contour of the sky component
gmm_oetscher_raster.png : the photograph with the patch grid drawn on it
gmm_oetscher_original.png : the photograph at patch resolution
gmm_oetscher_vegetation.png : brightness = posterior of the vegetation cluster
gmm_oetscher_mountain.png : brightness = posterior of the mountain cluster
gmm_oetscher_sky.png : brightness = posterior of the sky cluster
gmm_oetscher_vegetation_kmeans.png : the k-means assignment to the
vegetation cluster (bright = assigned)
gmm_oetscher_vegetation_sigma1.png, _sigma2.png, _sigma3.png :
brightness = posterior of the vegetation cluster with
every covariance matrix fixed at a small, a moderate
and a large multiple of the identity matrix
gmm.png : preview (checking only) -- the patches with the
three component ellipses drawn over them
"""
# /// script
# requires-python = ">=3.10"
# dependencies = ["matplotlib==3.10.7", "numpy==2.3.4"]
# ///
from pathlib import Path
import urllib.request
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.image import imread
OUT_DIR = Path(__file__).parent
ASSET_URL = "https://dictionaryofml.org/demoassets/"
def demo_asset(name):
"""Path to the input file `name`, fetched once if it is not already here.
Inside the repository the committed copy is used (pythondemos/ first,
then assets/), so a run there reproduces the committed outputs exactly.
Run from anywhere else -- from the copy published on dictionaryofml.org,
say -- the file is downloaded and cached beside the script, so the demo
is runnable without the repository.
"""
for cand in (OUT_DIR / name, OUT_DIR.parent / "assets" / name):
if cand.exists():
return cand
cached = OUT_DIR / name
print(f" fetching {name} from {ASSET_URL}{name}")
with urllib.request.urlopen(ASSET_URL + name, timeout=180) as resp:
cached.write_bytes(resp.read())
return cached
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}")
photograph cut into 27 x 36 = 972 patches of 128 x 128 pixels
Cut the photograph into 64x64 patches and measure how green and how blue each one is; check the patch count.
PATCH = 128
photo = imread(demo_asset("oetscher.jpg")) / 255.0
rows = photo.shape[0] // PATCH
cols = photo.shape[1] // PATCH
tiles = photo[:rows * PATCH, :cols * PATCH].reshape(
rows, PATCH, cols, PATCH, 3).mean(axis=(1, 3))
rgb = tiles.reshape(-1, 3)
green = rgb[:, 1] # average greenness, on a 0 to 1 scale
blue = rgb[:, 2] # average blueness, on the same scale
X = np.column_stack([green, blue])
with open(OUT_DIR / "gmm_patches.csv", "w") as f:
f.write("green,blue\n")
for a, b in X:
f.write(f"{a:.4f},{b:.4f}\n")
print(f" photograph cut into {rows} x {cols} = {len(X)} patches "
f"of {PATCH} x {PATCH} pixels")
check("[B-patches] every patch has a two-number feature vector",
X.shape == (rows * cols, 2))
[ok] [B-patches] every patch has a two-number feature vector
Fit the three-component GMM by the EM algorithm: E-step (posterior probability of each cluster index per patch), M-step (re-weighted cluster probabilities, means and covariance matrices). Check that the cluster probabilities sum to one and that the fit stops changing.
NRCLUSTER = 3
def normal_pdf(A, mean, cov):
"""Density of the multivariate normal distribution at each row of A."""
d = A.shape[1]
diff = A - mean
quad = np.einsum("ij,jk,ik->i", diff, np.linalg.inv(cov), diff)
return np.exp(-0.5 * quad) / np.sqrt(((2.0 * np.pi) ** d)
* np.linalg.det(cov))
def fit_quality(A, p, means, covs):
"""Negative log of the likelihood of the whole dataset."""
mix = sum(p[c] * normal_pdf(A, means[c], covs[c]) for c in range(len(p)))
return float(-np.log(mix).sum())
groups = np.array_split(np.argsort(X[:, 1]), NRCLUSTER) # ordered by blueness
means = np.stack([X[g].mean(axis=0) for g in groups])
covs = np.stack([np.cov(X[g].T) for g in groups])
p = np.full(NRCLUSTER, 1.0 / NRCLUSTER)
trace = [fit_quality(X, p, means, covs)]
for _ in range(300):
# E-step: posterior probability of each cluster index per patch
joint = np.stack([p[c] * normal_pdf(X, means[c], covs[c])
for c in range(NRCLUSTER)])
posterior = joint / joint.sum(axis=0)
# M-step: re-weighted cluster probabilities, means, covariance matrices
weight = posterior.sum(axis=1)
p = weight / len(X)
means = (posterior @ X) / weight[:, None]
covs = np.stack([
(posterior[c][:, None] * (X - means[c])).T @ (X - means[c]) / weight[c]
for c in range(NRCLUSTER)])
trace.append(fit_quality(X, p, means, covs))
order = np.argsort(means[:, 1]) # greenest first, sky last
p, means, covs = p[order], means[order], covs[order]
posterior = posterior[order]
label = posterior.argmax(axis=0)
check("[B-em] the cluster probabilities sum to one",
abs(p.sum() - 1.0) < 1e-12)
check("[B-em] the fit stops changing (last update tiny)",
abs(trace[-1] - trace[-2]) < 1e-9)
check("[B-em] no iteration increases the negative log-likelihood",
bool(np.all(np.diff(trace) <= 1e-9)))
[ok] [B-em] the cluster probabilities sum to one [ok] [B-em] the fit stops changing (last update tiny) [ok] [B-em] no iteration increases the negative log-likelihood vegetation probability 0.46, greenness +0.305, blueness +0.123, average color RGB [0.24 0.29 0.12], 15% of its patches in the upper half mountain and haze probability 0.39, greenness +0.401, blueness +0.496, average color RGB [0.26 0.43 0.54], 84% of its patches in the upper half bright sky probability 0.15, greenness +0.683, blueness +0.974, average color RGB [0.36 0.68 0.97], 100% of its patches in the upper half
The three components are the sunlit meadow, the dark foliage, and the sky with the mountain. Check each against the average color of its patches and against where they sit in the photograph, and check that each covariance matrix is symmetric positive semi-definite.
NAMES = ("vegetation", "mountain and haze", "bright sky")
where = np.argwhere(np.ones((rows, cols), dtype=bool)) # (row, col) per patch
for c in range(NRCLUSTER):
sel = label == c
print(f" {NAMES[c]:<14} probability {p[c]:.2f}, "
f"greenness {means[c, 0]:+.3f}, blueness {means[c, 1]:+.3f}, "
f"average color RGB {np.round(rgb[sel].mean(axis=0), 2)}, "
f"{(where[sel][:, 0] < rows / 2).mean():.0%} of its patches in the "
f"upper half")
check("[B-components] every average lies on the 0 to 1 scale",
float(X.min()) >= 0.0 and float(X.max()) <= 1.0)
check("[B-components] the sky is the bluest component",
means[2, 1] == means[:, 1].max())
check("[B-components] the sky sits in the upper half of the photograph",
float((where[label == 2][:, 0] < rows / 2).mean()) > 0.8)
check("[B-components] the vegetation is the least blue and sits in the "
"lower half", means[0, 1] == means[:, 1].min()
and float((where[label == 0][:, 0] >= rows / 2).mean()) > 0.8)
check("[B-components] each covariance matrix is symmetric",
all(np.allclose(c, c.T) for c in covs))
check("[B-components] each covariance matrix is positive semi-definite",
all(np.linalg.eigvalsh(c).min() > 0 for c in covs))
[ok] [B-components] every average lies on the 0 to 1 scale [ok] [B-components] the sky is the bluest component [ok] [B-components] the sky sits in the upper half of the photograph [ok] [B-components] the vegetation is the least blue and sits in the lower half [ok] [B-components] each covariance matrix is symmetric [ok] [B-components] each covariance matrix is positive semi-definite 820 of 972 patches are graded above 0.8 for one cluster; 152 are shared
Soft clustering: the posterior distribution grades membership. Check that most patches are graded decisively and that the rest are shared between clusters.
top = posterior.max(axis=0)
shared = top < 0.8
print(f" {int((~shared).sum())} of {len(X)} patches are graded above 0.8 "
f"for one cluster; {int(shared.sum())} are shared")
check("[B-soft] the posterior distribution of each patch sums to one",
np.allclose(posterior.sum(axis=0), 1.0))
check("[B-soft] most patches are graded decisively", (~shared).mean() > 0.8)
check("[B-soft] a shared patch is not assigned to exactly one cluster",
bool(np.all(posterior[:, shared].max(axis=0) < 0.8)))
[ok] [B-soft] the posterior distribution of each patch sums to one [ok] [B-soft] most patches are graded decisively [ok] [B-soft] a shared patch is not assigned to exactly one cluster
The photograph with the patch grid drawn on it, the photograph at patch resolution, and one copy per cluster whose brightness is that cluster's posterior probability. Because the posterior probabilities of a patch sum to one, the three copies add back up to the photograph.
def shrink_to_width(img, width_px):
"""Area-average an image down to width_px, keeping its aspect ratio.
Each output pixel is the mean of the input it covers, which is the right
filter for a photograph and needs nothing beyond numpy.
"""
if img.shape[1] <= width_px:
return np.asarray(img, dtype=float)
height_px = round(img.shape[0] * width_px / img.shape[1])
out = np.asarray(img, dtype=float)
for axis, n_out in ((1, width_px), (0, height_px)):
n_in = out.shape[axis]
edges = np.linspace(0.0, n_in, n_out + 1)
cum = np.cumsum(out, axis=axis)
cum = np.concatenate(
[np.zeros_like(np.take(cum, [0], axis=axis)), cum], axis=axis)
grid = np.arange(n_in + 1)
cum = np.apply_along_axis(
lambda col, g=grid, e=edges: np.interp(e, g, col), axis, cum)
shape = [1] * out.ndim
shape[axis] = n_out
out = np.diff(cum, axis=axis) / np.diff(edges).reshape(shape)
return out
def save_image(arr, name, zoom=6):
"""Write an RGB array as a PNG, enlarged so the patches stay visible."""
img = np.clip(arr, 0.0, 1.0).repeat(zoom, axis=0).repeat(zoom, axis=1)
plt.imsave(OUT_DIR / name, img)
# the photograph at full detail with the patch grid drawn on it, so
# the reader can see what one data point covers
# The raster is printed at 0.56\textwidth, and the book's textwidth is 5.40 in,
# so it is 3.02 in across on the page; at 300 ppi that is 907 px against the
# 1152 px of the full-detail view, which went into the PDF at over 400 ppi.
# The photograph is shrunk first and the patch rules drawn afterwards, so every
# rule stays exactly one white pixel; shrinking the ruled image would blur the
# rules, which are what shows the reader the extent of one data point.
SHRINK, TARGET_DPI = 4, 300
RASTER_WIDTH_IN = 3.02
view = photo[:rows * PATCH, :cols * PATCH:, :][::SHRINK, ::SHRINK]
cell = -(-round(RASTER_WIDTH_IN * TARGET_DPI) // cols) # ceil: at least 300 ppi
view = shrink_to_width(view, cols * cell)
view[::cell, :, :] = 1.0 # horizontal rules
view[:, ::cell, :] = 1.0 # vertical rules
view[-1, :, :] = 1.0
view[:, -1, :] = 1.0
save_image(view, "gmm_oetscher_raster.png", zoom=1)
check("[B-images] the raster has one cell per patch",
view.shape[:2] == (rows * cell, cols * cell))
save_image(tiles, "gmm_oetscher_original.png")
COPIES = ("gmm_oetscher_vegetation.png", "gmm_oetscher_mountain.png",
"gmm_oetscher_sky.png")
dimmed = [tiles * posterior[c].reshape(rows, cols, 1)
for c in range(NRCLUSTER)]
for fname, arr in zip(COPIES, dimmed):
save_image(arr, fname)
check("[B-images] one copy of the photograph per cluster",
all((OUT_DIR / f).exists() for f in COPIES))
check("[B-images] the three copies add back up to the photograph",
np.allclose(sum(dimmed), tiles))
[ok] [B-images] the raster has one cell per patch [ok] [B-images] one copy of the photograph per cluster [ok] [B-images] the three copies add back up to the photograph variance 0.01: 85% of posteriors above 0.99, assignment agrees with k-means on 99.8% of patches variance 1e-05: 100% of posteriors above 0.99, assignment agrees with k-means on 99.8% of patches
With equal cluster probabilities and covariance matrices that shrink to a multiple of the identity matrix, the posterior concentrates on the nearest mean and the EM update becomes the k-means update, so the algorithm reduces to k-means. Check both.
def kmeans(A, cents):
"""Lloyd's algorithm from the given starting means."""
while True:
lab = ((A[:, None, :] - cents) ** 2).sum(axis=2).argmin(axis=1)
new = np.stack([A[lab == c].mean(axis=0) for c in range(len(cents))])
if np.allclose(new, cents):
return lab, new
cents = new
def spherical_em(A, cents, var):
"""EM with equal cluster probabilities and covariance matrices var*I."""
post = None
for _ in range(300):
d2 = ((A[:, None, :] - cents) ** 2).sum(axis=2)
logp = -d2 / (2.0 * var)
post = np.exp(logp - logp.max(axis=1, keepdims=True))
post /= post.sum(axis=1, keepdims=True)
new = (post.T @ A) / post.sum(axis=0)[:, None]
if np.allclose(new, cents):
break
cents = new
return post.argmax(axis=1), cents, post
start = np.stack([X[g].mean(axis=0) for g in groups])
lab_km, cent_km = kmeans(X, start.copy())
for var in (1e-2, 1e-5):
lab_em, cent_em, post = spherical_em(X, start.copy(), var)
agree = float((lab_em == lab_km).mean())
hard = float((post.max(axis=1) > 0.99).mean())
print(f" variance {var:g}: {hard:.0%} of posteriors above 0.99, "
f"assignment agrees with k-means on {agree:.1%} of patches")
if var == 1e-5:
check("[B-kmeans] the posterior becomes a hard assignment", hard > 0.99)
check("[B-kmeans] the assignment agrees with k-means", agree > 0.99)
check("[B-kmeans] the means agree with the k-means centroids",
np.abs(cent_em - cent_km).max() < 0.01)
[ok] [B-kmeans] the posterior becomes a hard assignment [ok] [B-kmeans] the assignment agrees with k-means [ok] [B-kmeans] the means agree with the k-means centroids variance 0.0001: vegetation posterior between 0.00 and 1.00, 0% of patches graded between 0.2 and 0.8 variance 0.02: vegetation posterior between 0.00 and 0.78, 63% of patches graded between 0.2 and 0.8 variance 1: vegetation posterior between 0.33 and 0.33, 100% of patches graded between 0.2 and 0.8
The same constrained model for a small, a moderate and a large common variance: the photograph dimmed by the posterior probability of the vegetation cluster, next to the k-means assignment to that cluster. Check that the small variance reproduces the k-means assignment, the moderate one grades the patches, and the large one grades every patch equally.
SIGMAS = (1e-4, 2e-2, 1.0) # small, moderate, large common variance
mask_km = (lab_km == 0).astype(float).reshape(rows, cols, 1)
save_image(tiles * mask_km, "gmm_oetscher_vegetation_kmeans.png")
graded = []
for i, var in enumerate(SIGMAS, start=1):
_, _, post = spherical_em(X, start.copy(), var)
veg = post[:, 0] # posterior of the vegetation cluster
graded.append(veg)
save_image(tiles * veg.reshape(rows, cols, 1),
f"gmm_oetscher_vegetation_sigma{i}.png")
print(f" variance {var:g}: vegetation posterior between {veg.min():.2f} "
f"and {veg.max():.2f}, {((veg > 0.2) & (veg < 0.8)).mean():.0%} of "
f"patches graded between 0.2 and 0.8")
check("[B-sigma] the small variance reproduces the k-means assignment",
float(((graded[0] > 0.5) == (lab_km == 0)).mean()) > 0.99)
check("[B-sigma] the moderate variance grades a sizeable share of patches",
float(((graded[1] > 0.2) & (graded[1] < 0.8)).mean()) > 0.3)
check("[B-sigma] the large variance grades every patch equally",
np.allclose(graded[2], 1.0 / NRCLUSTER, atol=1e-3))
[ok] [B-sigma] the small variance reproduces the k-means assignment [ok] [B-sigma] the moderate variance grades a sizeable share of patches [ok] [B-sigma] the large variance grades every patch equally
Write the patch scatter and the three component ellipses for the entry's figure, plus the preview.
def ellipse(mean, cov, sigma=2.0, n=200):
"""Contour at `sigma` standard deviations of one component."""
vals, vecs = np.linalg.eigh(cov)
t = np.linspace(0.0, 2.0 * np.pi, n)
circle = np.stack([np.cos(t), np.sin(t)])
return (mean[:, None] + sigma * vecs @ (np.sqrt(vals)[:, None]
* circle)).T
with open(OUT_DIR / "gmm_points.csv", "w") as f:
f.write("x1,x2\n")
for a, b in X:
f.write(f"{a:.4f},{b:.4f}\n")
for c, name in enumerate(("gmm_ellipse1.csv", "gmm_ellipse2.csv",
"gmm_ellipse3.csv")):
with open(OUT_DIR / name, "w") as f:
f.write("x1,x2\n")
for a, b in ellipse(means[c], covs[c]):
f.write(f"{a:.4f},{b:.4f}\n")
fig, ax = plt.subplots(figsize=(6.2, 5.0))
ax.plot(X[:, 0], X[:, 1], "o", color="0.6", markersize=2.2,
linestyle="none", label="patch")
for c, style in enumerate(("-", "--", ":")):
e = ellipse(means[c], covs[c])
ax.plot(e[:, 0], e[:, 1], color="black", linewidth=1.7, linestyle=style,
label=NAMES[c])
ax.set_xlabel("average greenness of the patch")
ax.set_ylabel("average blueness of the patch")
ax.set_title("Patches of the Oetscher photograph and a three-component GMM")
ax.legend(frameon=False, loc="upper right", fontsize="small")
fig.tight_layout()
fig.savefig(OUT_DIR / "gmm.png", dpi=150)
plt.close(fig)
check("[B-plot] the three ellipses and the patch scatter were written",
(OUT_DIR / "gmm_ellipse3.csv").exists()
and (OUT_DIR / "gmm_points.csv").exists())
passed = sum(1 for _, ok in report if ok)
print(f"\n{passed}/{len(report)} checks pass")
[ok] [B-plot] the three ellipses and the patch scatter were written 23/23 checks pass
B-plot writes when the script runs