P-picture writes when the script runsDictionary of Applied Machine Learning · probability distribution
Numerical companion to the entry probability distribution: it recomputes what the entry states and prints one line per check
One block per paragraph of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts. Self-contained (numpy/matplotlib only), fixed seed.
Run it without installing anything:uv run https://dictionaryofml.org/terms/probdist.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 probdist.py into a folder and run uv run probdist.py there. With NumPy and Matplotlib already installed, python3 probdist.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download probdist.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""
probdist.py — numerical companion to the glossary entry
'probability distribution'.
One block per paragraph of the entry (marked [P...]): each block verifies
numerically what the corresponding statement asserts. Self-contained
(numpy/matplotlib only), fixed seed.
The method used throughout the first blocks is the simplest one
available: it reads a data set of numbers and delivers their average.
That keeps every quantity below available in closed form, so a check
compares a measured number against a formula rather than against
another simulation.
Blocks
------
[P-method] One distribution, many data sets, a different output each
time: the outputs of the method are realizations of another
RV, whose spread shrinks as sigma/sqrt(m) with the data set
size m. The distribution of that output RV depends on the
common distribution P and on the method: a different method
(taking the middle value of the sorted data set) run on the
same data sets has a larger output spread.
[P-typical] The distribution decides which data points are typical
(relative frequencies of a growing sample converge to it),
and it can be estimated from a data set: the empirical
frequencies computed from one large data set recover the
distribution that generated it.
[P-specify] How a distribution is specified: a binary RV by the single
probability P(y = 0), and a continuous real-valued RV by a
pdf p, for which P(x in [a, b]) ~ p(a)|b - a| on a short
interval.
[P-picture] The distribution visualized: the joint pdf of a data point
(feature, label) as a grayscale over the plane, a
two-component Gaussian mixture. Paralleling the entry's
Fig. 1, three data sets of realizations of iid RVs with
this distribution are drawn, and three hypothesis maps are
learned from them by the same polynomial regression method.
The checks verify that the plotted window carries the
probability mass, that the data sets concentrate where the
pdf is large, that the three learned maps differ, and that
they differ most where the pdf is small.
Outputs
-------
probdist.png : preview figure (checking only).
probdist_density.csv : joint pdf on a grid (columns x, y, p; row-wise
in y).
probdist_train1.csv, probdist_train2.csv, probdist_train3.csv
: the three data sets (columns x, y).
probdist_hypotheses.csv : the three learned hypothesis curves
(columns x, yhat1, yhat2, yhat3).
Data generated by pythondemos/probdist.py.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from pathlib import Path
OUT_DIR = Path(__file__).parent
rng = np.random.default_rng(42)
report = []
MU, SIGMA = 1.0, 2.0 # the common distribution: N(MU, SIGMA^2)
NR_DATASETS = 4000
def check(name, ok):
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
def method(data):
"""The ML method: it delivers the average of the data set."""
return float(np.mean(data))
One distribution, many data sets, a different output each time: the outputs of the method are realizations of another RV, whose spread shrinks as sigma/sqrt(m) with the data set size m. The distribution of that output RV depends on the common distribution P and on the method: a different method (taking the middle value of the sorted data set) run on the same data sets has a larger output spread.
print("[P-method] one distribution, many data sets, a different output each time")
outputs = {}
for m in (25, 400):
outputs[m] = np.array([method(rng.normal(MU, SIGMA, m))
for _ in range(NR_DATASETS)])
spread = {m: float(np.std(out)) for m, out in outputs.items()}
print(f" output spread: {spread[25]:.4f} at m=25, {spread[400]:.4f} at m=400"
f" (sigma/sqrt(m) = {SIGMA / np.sqrt(25):.4f}, "
f"{SIGMA / np.sqrt(400):.4f})")
check("two data sets from the same distribution give different outputs",
outputs[25][0] != outputs[25][1])
check("the spread of the outputs shrinks with the data set size",
spread[400] < spread[25])
for m in (25, 400):
check(f"the spread matches sigma/sqrt(m) within 5% at m={m}",
abs(spread[m] - SIGMA / np.sqrt(m)) / (SIGMA / np.sqrt(m)) < 0.05)
# the distribution of the output RV depends on P and on the method: the
# middle value of a sorted Gaussian data set has a larger spread than
# the average
med_outputs = np.array([float(np.median(rng.normal(MU, SIGMA, 25)))
for _ in range(NR_DATASETS)])
print(f" output spread at m=25: {spread[25]:.4f} (average), "
f"{float(np.std(med_outputs)):.4f} (middle value)")
check("a different method (the middle value) has a different output "
"distribution: its spread is larger",
float(np.std(med_outputs)) > 1.1 * spread[25])
[P-method] one distribution, many data sets, a different output each time
output spread: 0.4034 at m=25, 0.0979 at m=400 (sigma/sqrt(m) = 0.4000, 0.1000)
[ok] two data sets from the same distribution give different outputs
[ok] the spread of the outputs shrinks with the data set size
[ok] the spread matches sigma/sqrt(m) within 5% at m=25
[ok] the spread matches sigma/sqrt(m) within 5% at m=400
output spread at m=25: 0.4034 (average), 0.4911 (middle value)
[ok] a different method (the middle value) has a different output distribution: its spread is larger
The distribution decides which data points are typical (relative frequencies of a growing sample converge to it), and it can be estimated from a data set: the empirical frequencies computed from one large data set recover the distribution that generated it.
print("[P-typical] typical data points; the distribution estimated from a data set")
p_true = np.array([0.5, 0.3, 0.2]) # distribution on {0, 1, 2}
errs = []
for m in (10**2, 10**4, 10**6):
draws = rng.choice(3, size=m, p=p_true)
errs.append(np.max(np.abs(np.bincount(draws, minlength=3) / m - p_true)))
print(f" max |frequency - p| for m=1e2,1e4,1e6: "
f"{errs[0]:.4f}, {errs[1]:.4f}, {errs[2]:.4f}")
check("relative frequencies converge to the distribution", errs[0] > errs[2])
check("the empirical frequencies of one large data set estimate the "
"distribution that generated it (within 2e-3)", errs[2] < 2e-3)
[P-typical] typical data points; the distribution estimated from a data set
max |frequency - p| for m=1e2,1e4,1e6: 0.0300, 0.0110, 0.0005
[ok] relative frequencies converge to the distribution
[ok] the empirical frequencies of one large data set estimate the distribution that generated it (within 2e-3)
How a distribution is specified: a binary RV by the single probability P(y = 0), and a continuous real-valued RV by a pdf p, for which P(x in [a, b]) ~ p(a)|b - a| on a short interval.
print("[P-specify] one probability specifies a binary RV; a pdf a continuous one")
p0 = 0.73
y = (rng.uniform(size=10**6) >= p0).astype(int) # P(y = 0) = p0
f0 = float(np.mean(y == 0))
check("empirical P(y = 0) recovers p0 = 0.73", abs(f0 - p0) < 2e-3)
check("P(y = 1) = 1 - P(y = 0)", np.isclose(np.mean(y == 1), 1 - f0))
pdf = lambda t: np.exp(-t ** 2 / 2) / np.sqrt(2 * np.pi)
x = rng.standard_normal(10**7)
a = 0.5
rel_errs = []
for width in (0.5, 0.1, 0.02):
p_emp = np.mean((x >= a) & (x <= a + width))
rel_errs.append(abs(p_emp - pdf(a) * width) / p_emp)
print(f" relative approximation error for |b-a|=0.5,0.1,0.02: "
f"{rel_errs[0]:.3f}, {rel_errs[1]:.3f}, {rel_errs[2]:.3f}")
check("the pdf approximation improves as the interval shrinks",
rel_errs[0] > rel_errs[1] > rel_errs[2])
check("relative error below 1% for |b - a| = 0.02", rel_errs[2] < 0.01)
[P-specify] one probability specifies a binary RV; a pdf a continuous one
[ok] empirical P(y = 0) recovers p0 = 0.73
[ok] P(y = 1) = 1 - P(y = 0)
relative approximation error for |b-a|=0.5,0.1,0.02: 0.174, 0.027, 0.007
[ok] the pdf approximation improves as the interval shrinks
[ok] relative error below 1% for |b - a| = 0.02
The distribution visualized: the joint pdf of a data point (feature, label) as a grayscale over the plane, a two-component Gaussian mixture. Paralleling the entry's Fig. 1, three data sets of realizations of iid RVs with this distribution are drawn, and three hypothesis maps are learned from them by the same polynomial regression method. The checks verify that the plotted window carries the probability mass, that the data sets concentrate where the pdf is large, that the three learned maps differ, and that they differ most where the pdf is small.
print("[P-picture] grayscale pdf; three data sets, three learned maps")
# The joint pdf of a data point z = (x, y): the feature x follows a
# two-component Gaussian mixture, and the label y is the value of a fixed
# curve at x plus Gaussian noise.
MIX_W = np.array([0.5, 0.5])
MIX_MU = np.array([-2.0, 2.0])
MIX_S = np.array([0.7, 0.9])
NOISE = 0.35
M_TRAIN, NR_SETS = 40, 3
def gauss_pdf(t, mu, s):
return np.exp(-((t - mu) ** 2) / (2 * s ** 2)) / (s * np.sqrt(2 * np.pi))
def pdf_feature(t):
return sum(w * gauss_pdf(t, mu, s) for w, mu, s in zip(MIX_W, MIX_MU, MIX_S))
def pdf_joint(t, u):
return pdf_feature(t) * gauss_pdf(u, np.tanh(t), NOISE)
def draw_points(m):
comp = rng.choice(2, size=m, p=MIX_W)
xx = rng.normal(MIX_MU[comp], MIX_S[comp])
yy = np.tanh(xx) + NOISE * rng.standard_normal(m)
return xx, yy
# three data sets of iid draws, and the map the same polynomial regression
# method learns from each (the degree-3 polynomial minimizing the average
# squared deviation on the data set, available in closed form)
x_curve = np.linspace(-4.4, 4.4, 200)
sets, curves = [], []
for _ in range(NR_SETS):
x_tr, y_tr = draw_points(M_TRAIN)
design = np.vander(x_tr, 4)
w_hat = np.linalg.solve(design.T @ design, design.T @ y_tr)
sets.append((x_tr, y_tr))
curves.append(np.polyval(w_hat, x_curve))
curves = np.array(curves)
xs = np.linspace(-4.4, 4.4, 61)
ys = np.linspace(-1.9, 1.9, 41)
grid_x, grid_y = np.meshgrid(xs, ys)
grid_p = pdf_joint(grid_x, grid_y)
mass = float(grid_p.sum() * (xs[1] - xs[0]) * (ys[1] - ys[0]))
print(f" probability mass inside the plotted window: {mass:.4f}")
check("the plotted window carries the probability mass (within 2%)",
abs(mass - 1.0) < 0.02)
for i, (x_tr, y_tr) in enumerate(sets, 1):
check(f"data set {i} concentrates where the pdf is large",
float(pdf_joint(x_tr, y_tr).mean()) > float(grid_p.mean()))
gap01 = float(np.max(np.abs(curves[0] - curves[1])))
check("the three learned maps differ", gap01 > 0.05)
spread_curve = curves.std(axis=0) # pointwise over the maps
dens = pdf_feature(x_curve)
low, high = dens < np.median(dens), dens >= np.median(dens)
print(f" map spread: {spread_curve[low].mean():.3f} where the pdf is "
f"small, {spread_curve[high].mean():.3f} where it is large")
check("the maps differ most where the pdf is small",
spread_curve[low].mean() > spread_curve[high].mean())
np.savetxt(OUT_DIR / "probdist_density.csv",
np.column_stack([grid_x.ravel(), grid_y.ravel(), grid_p.ravel()]),
delimiter=",", header="x,y,p", comments="", fmt="%.6f")
for i, (x_tr, y_tr) in enumerate(sets, 1):
np.savetxt(OUT_DIR / f"probdist_train{i}.csv",
np.column_stack([x_tr, y_tr]),
delimiter=",", header="x,y", comments="", fmt="%.6f")
np.savetxt(OUT_DIR / "probdist_hypotheses.csv",
np.column_stack([x_curve, curves[0], curves[1], curves[2]]),
delimiter=",", header="x,yhat1,yhat2,yhat3", comments="",
fmt="%.6f")
# ------------------------------------------------------------ preview
fig, ax = plt.subplots(1, 3, figsize=(12.6, 3.0))
for m, style in ((25, "--"), (400, "-")):
ax[0].hist(outputs[m], bins=60, histtype="step", density=True,
color="k", linestyle=style, label=f"m = {m}")
ax[0].set_xlabel("output of the method (the average)")
ax[0].set_ylabel("density over data sets")
ax[0].set_title("outputs of one method over 4000 data sets")
ax[0].legend(frameon=False)
t = np.linspace(-4, 4, 400)
ax[1].plot(t, pdf(t), "k-")
ax[1].fill_between(t, pdf(t), where=(t >= a) & (t <= a + 0.5),
facecolor="none", hatch="///", edgecolor="k")
ax[1].set_xlabel("value of the RV")
ax[1].set_ylabel("probability density p")
ax[1].set_title("P(x in [a, b]) is the shaded area")
ax[2].pcolormesh(xs, ys, grid_p, cmap="Greys", shading="nearest",
vmin=0.0, vmax=float(grid_p.max()) * 1.3)
for i, (ls, mk) in enumerate(zip(("-", "--", ":"), ("o", "s", "^"))):
ax[2].plot(x_curve, curves[i], "k", linestyle=ls,
label=f"learned map {i + 1}")
ax[2].scatter(sets[i][0], sets[i][1], s=18, marker=mk,
facecolors=("black", "white", "0.5")[i],
edgecolors=("white", "black", "black")[i],
linewidths=0.4, label=f"data set {i + 1}", zorder=3)
ax[2].set_xlabel("feature x")
ax[2].set_ylabel("label y")
ax[2].set_title("pdf (grayscale), three data sets, three learned maps")
ax[2].legend(frameon=False, loc="upper left", fontsize=7)
fig.tight_layout()
fig.savefig(OUT_DIR / "probdist.png", dpi=110)
print(f"\n{sum(ok for _, ok in report)}/{len(report)} checks passed")
assert all(ok for _, ok in report)
[P-picture] grayscale pdf; three data sets, three learned maps
probability mass inside the plotted window: 0.9963
[ok] the plotted window carries the probability mass (within 2%)
[ok] data set 1 concentrates where the pdf is large
[ok] data set 2 concentrates where the pdf is large
[ok] data set 3 concentrates where the pdf is large
[ok] the three learned maps differ
map spread: 0.105 where the pdf is small, 0.059 where it is large
[ok] the maps differ most where the pdf is small
17/17 checks passed
P-picture writes when the script runs