diff --git a/docs/source/_gen_diagrams.py b/docs/source/_gen_diagrams.py
new file mode 100644
index 0000000..07ed749
--- /dev/null
+++ b/docs/source/_gen_diagrams.py
@@ -0,0 +1,688 @@
+#!/usr/bin/env python3
+"""
+Generate all matplotlib diagrams for mathematical_foundations.md.
+
+Run from the docs/source directory (or workspace root):
+ python docs/source/_gen_diagrams.py
+
+Outputs SVG files to docs/source/_static/diagrams/
+"""
+
+import os
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import matplotlib.ticker as mticker
+from scipy.stats import norm
+
+# ─── output dir ─────────────────────────────────────────────────────────────
+OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "_static", "diagrams")
+os.makedirs(OUT, exist_ok=True)
+
+# ─── palette & defaults ─────────────────────────────────────────────────────
+C0 = "#2E6BE5" # blue
+C1 = "#E8850A" # orange
+C2 = "#27AE60" # green
+C3 = "#D62728" # red
+GRAY = "#888888"
+BAND = "#AACBE8"
+
+matplotlib.rcParams.update({
+ "font.size" : 11,
+ "axes.titlesize" : 12,
+ "axes.labelsize" : 11,
+ "xtick.labelsize" : 9,
+ "ytick.labelsize" : 9,
+ "axes.spines.top" : False,
+ "axes.spines.right" : False,
+ "figure.dpi" : 150,
+ "savefig.bbox" : "tight",
+ "savefig.transparent" : False,
+ "figure.facecolor" : "white",
+ "axes.facecolor" : "white",
+ "lines.linewidth" : 1.8,
+ "text.usetex" : False,
+})
+
+def save(name):
+ plt.savefig(os.path.join(OUT, name + ".svg"))
+ plt.close()
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §1 DIFFERENTIAL EVOLUTION
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_de_mutation():
+ r1 = np.array([0.5, 0.3])
+ r2 = np.array([1.2, 1.4])
+ r3 = np.array([1.8, 0.6])
+ F = 0.7
+ vi = r1 + F * (r2 - r3)
+
+ fig, ax = plt.subplots(figsize=(6, 4.2))
+
+ # difference vector r3 → r2
+ ax.annotate("", r2, r3,
+ arrowprops=dict(arrowstyle="-|>", color=C2, lw=2.0, mutation_scale=14))
+ mid = (r2 + r3) / 2
+ ax.text(mid[0] - 0.05, mid[1] + 0.09,
+ r"$F(\mathbf{x}_{r_2}-\mathbf{x}_{r_3})$",
+ ha="center", fontsize=10, color=C2)
+
+ # mutation arrow r1 → vi (dashed)
+ ax.annotate("", vi, r1,
+ arrowprops=dict(arrowstyle="-|>", color=C1, lw=2.0,
+ mutation_scale=14, linestyle="dashed"))
+ ax.text((r1[0]+vi[0])/2, (r1[1]+vi[1])/2 - 0.1,
+ r"$+F(\cdots)$", ha="center", fontsize=9, color=C1)
+
+ pts = {
+ r"$\mathbf{x}_{r_1}$ (base)": (r1, C0),
+ r"$\mathbf{x}_{r_2}$": (r2, C0),
+ r"$\mathbf{x}_{r_3}$": (r3, C0),
+ r"$\mathbf{v}_i$ (mutant)": (vi, C1),
+ }
+ for lbl, (p, col) in pts.items():
+ ax.scatter(*p, s=90, color=col, zorder=6)
+ offset = (0.05, 0.07)
+ if "mutant" in lbl:
+ offset = (0.07, 0.05)
+ ax.text(p[0] + offset[0], p[1] + offset[1], lbl, fontsize=10, color=col)
+
+ ax.set_xlim(0.1, 2.5); ax.set_ylim(0.0, 1.85)
+ ax.set_xlabel(r"$x_1$"); ax.set_ylabel(r"$x_2$")
+ ax.set_title(r"DE Mutation: $\mathbf{v}_i = \mathbf{x}_{r_1} + F\,(\mathbf{x}_{r_2} - \mathbf{x}_{r_3})$")
+ ax.set_aspect("equal", adjustable="box")
+ save("fig_de_mutation")
+
+
+def fig_rastrigin():
+ x = np.linspace(-2.5, 2.5, 800)
+ y = 10 + x**2 - 10 * np.cos(2 * np.pi * x)
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(x, y, color=C0, lw=2, label=r"$f(x) = 10 + x^2 - 10\cos(2\pi x)$")
+ ax.fill_between(x, y, alpha=0.07, color=C0)
+ ax.axhline(0, color=GRAY, lw=0.7, ls=":")
+
+ # global minimum
+ ax.scatter([0], [0], s=110, color=C1, zorder=6, label=r"global min $f^*=0$", marker="*")
+
+ # local minima
+ lm_x = np.array([-2.0, -1.0, 1.0, 2.0])
+ lm_y = 10 + lm_x**2 - 10 * np.cos(2 * np.pi * lm_x)
+ ax.scatter(lm_x, lm_y, s=55, color=C3, zorder=5, label="local minima", marker="o")
+
+ ax.annotate(r"$\approx 10^d$ local pits", (1.0, lm_y[2]),
+ (1.5, 12), fontsize=9, color=C3,
+ arrowprops=dict(arrowstyle="->", color=C3, lw=1.0))
+
+ ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$f(x)$")
+ ax.set_title(r"Rastrigin function ($d = 1$) — many local minima")
+ ax.legend(fontsize=9, framealpha=0.6)
+ save("fig_rastrigin")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.1 BROWNIAN MOTION
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_random_walk():
+ rng = np.random.default_rng(42)
+ n = 300
+ t = np.linspace(0, 1, n)
+ W = np.cumsum(rng.choice([-1, 1], size=n)) / np.sqrt(n)
+
+ fig, ax = plt.subplots(figsize=(7, 3.5))
+ ax.plot(t, W, color=C0, lw=1.4)
+ ax.axhline(0, color=GRAY, lw=0.8, ls="--", alpha=0.6)
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$W_t^{(n)}$")
+ ax.set_title(r"Coin-flip random walk ($n=300$) $\longrightarrow$ Brownian motion as $n\to\infty$")
+ save("fig_random_walk")
+
+
+def fig_bm_fan():
+ rng = np.random.default_rng(0)
+ n, dt = 500, 0.002
+ npaths = 10
+ ts = np.linspace(0, 1, n)
+ paths = np.cumsum(rng.normal(0, np.sqrt(dt), (npaths, n)), axis=1)
+ paths[:, 0] = 0
+
+ fig, ax = plt.subplots(figsize=(7, 4.2))
+ lo, hi = -2 * np.sqrt(ts), 2 * np.sqrt(ts)
+ ax.fill_between(ts, lo, hi, alpha=0.13, color=C0, label=r"$\pm 2\sqrt{t}$ (95% band)")
+ ax.plot(ts, hi, color=C0, lw=1.2, ls="--", alpha=0.55)
+ ax.plot(ts, lo, color=C0, lw=1.2, ls="--", alpha=0.55)
+ colors_cycle = plt.cm.tab10(np.linspace(0, 0.9, npaths))
+ for i, p in enumerate(paths):
+ ax.plot(ts, p, lw=0.9, alpha=0.75, color=colors_cycle[i])
+ ax.axhline(0, color=GRAY, lw=0.8, ls=":")
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$W_t$")
+ ax.set_title(r"Brownian motion — sample paths spread as $\sqrt{t}$ (trumpet fan)")
+ ax.legend(fontsize=9, framealpha=0.7)
+ save("fig_bm_fan")
+
+
+def fig_gbm():
+ rng = np.random.default_rng(7)
+ T, n, dt = 1.0, 500, 0.002
+ mu, sigma, S0 = 0.10, 0.30, 1.0
+ ts = np.linspace(0, T, n)
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(ts, S0 * np.exp(mu * ts), color=C1, lw=1.8, ls="--",
+ label=r"$\mathbb{E}[S_t] = S_0 e^{\mu t}$")
+ ax.plot(ts, S0 * np.exp((mu - 0.5*sigma**2) * ts), color=C2, lw=1.5, ls=":",
+ label=r"median $\approx S_0 e^{(\mu-\sigma^2/2)t}$")
+ colors_cycle = plt.cm.Blues(np.linspace(0.4, 0.85, 7))
+ for i in range(7):
+ W = np.cumsum(rng.normal(0, np.sqrt(dt), n))
+ S = S0 * np.exp((mu - 0.5*sigma**2) * ts + sigma * W)
+ ax.plot(ts, S, lw=0.9, alpha=0.7, color=colors_cycle[i])
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$S_t$")
+ ax.set_title(r"Geometric Brownian motion ($\mu=0.10,\;\sigma=0.30$)")
+ ax.legend(fontsize=9, framealpha=0.6)
+ save("fig_gbm")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.2 ITŌ CALCULUS
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_ito_correction():
+ t = np.linspace(0, 2.2, 300)
+ mu, sigma = 0.12, 0.30
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(t, mu * t, color=C1, lw=2, ls="--",
+ label=r"Naïve slope $\mu t$ (wrong)")
+ ax.plot(t, (mu - 0.5*sigma**2) * t, color=C0, lw=2,
+ label=r"Itō slope $(\mu - \sigma^2/2)\,t$ (correct)")
+
+ # gap annotation at t = 1.8
+ g_x = 1.8
+ y_top = mu * g_x
+ y_bot = (mu - 0.5*sigma**2) * g_x
+ ax.annotate("", (g_x, y_bot), (g_x, y_top),
+ arrowprops=dict(arrowstyle="<->", color=C3, lw=1.6))
+ ax.text(g_x + 0.07, (y_top + y_bot) / 2,
+ r"gap $= \sigma^2 T/2$", fontsize=9, color=C3, va="center")
+
+ ax.axhline(0, color=GRAY, lw=0.6, ls=":")
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$\mathbb{E}[\log S_t] - \log S_0$")
+ ax.set_title(r"Itō correction: $\mathbb{E}[\log S_t]$ always below the naïve slope $\mu t$")
+ ax.legend(fontsize=9)
+ save("fig_ito_correction")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.3 FOKKER-PLANCK
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_fokker_planck():
+ x = np.linspace(-0.5, 5.5, 600)
+ mu_drift, sigma_diff = 0.8, 0.3
+ times = [0.05, 0.5, 1.5]
+ colors = [C3, C2, C0]
+ labels = [r"$t = 0.05$ (narrow spike)",
+ r"$t = 0.50$",
+ r"$t = 1.50$ (wide, drifted)"]
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ for t, col, lbl in zip(times, colors, labels):
+ mean = mu_drift * t
+ std = sigma_diff * np.sqrt(t)
+ y = norm.pdf(x, mean, std)
+ ax.plot(x, y, color=col, lw=2, label=lbl)
+ ax.fill_between(x, y, alpha=0.10, color=col)
+
+ ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$p(t, x)$")
+ ax.set_title(r"Fokker-Planck: density drifts $(\mu=0.8)$ and broadens $(\sigma=0.3)$")
+ ax.legend(fontsize=9)
+ save("fig_fokker_planck")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.3 EULER-MARUYAMA vs MILSTEIN
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_em_milstein():
+ dts = np.array([0.1, 0.05, 0.02, 0.01, 0.005, 0.001])
+ em_err = 0.38 * dts**0.5
+ mil_err = 0.19 * dts**1.0
+
+ fig, ax = plt.subplots(figsize=(6, 4))
+ ax.loglog(dts, em_err, "o-", color=C0, lw=2, ms=7,
+ label=r"Euler-Maruyama (order $1/2$)")
+ ax.loglog(dts, mil_err, "s--", color=C1, lw=2, ms=7,
+ label=r"Milstein (order $1$)")
+ ax.set_xlabel(r"Step size $\Delta t$")
+ ax.set_ylabel(r"Strong error $\|X_T - \hat{X}_T\|$")
+ ax.set_title("SDE numerical schemes — strong convergence order")
+ ax.legend(fontsize=10); ax.grid(True, which="both", alpha=0.3)
+ save("fig_em_milstein")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.4 ORNSTEIN-UHLENBECK
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_ou_path():
+ rng = np.random.default_rng(3)
+ T, n, dt = 5.0, 2000, 0.0025
+ kappa, theta, sigma = 3.0, 0.5, 0.4
+ X = np.zeros(n); X[0] = 2.0
+ for i in range(1, n):
+ X[i] = X[i-1] + kappa * (theta - X[i-1]) * dt + sigma * rng.normal(0, np.sqrt(dt))
+
+ ts = np.linspace(0, T, n)
+ sig_inf = sigma / np.sqrt(2 * kappa)
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(ts, X, color=C0, lw=1.0, alpha=0.9, label=r"$X_t$")
+ ax.axhline(theta, color=C1, lw=1.8, ls="--",
+ label=fr"$\theta = {theta}$ (long-run mean)")
+ ax.fill_between(ts,
+ theta - 2 * sig_inf,
+ theta + 2 * sig_inf,
+ alpha=0.10, color=GRAY, label=r"$\theta \pm 2\sigma_\infty$")
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$X_t$")
+ ax.set_title(fr"Ornstein-Uhlenbeck ($\kappa={kappa},\;\theta={theta},\;\sigma={sigma}$) — mean-reversion")
+ ax.legend(fontsize=9)
+ save("fig_ou_path")
+
+
+def fig_ou_transition():
+ x = np.linspace(-0.3, 2.6, 500)
+ kappa, theta, sigma, x0 = 3.0, 0.5, 0.4, 2.0
+ taus = [0.1, 0.5, 2.0]
+ colors = [C3, C2, C0]
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ for tau, col in zip(taus, colors):
+ mean = theta + (x0 - theta) * np.exp(-kappa * tau)
+ var = sigma**2 / (2 * kappa) * (1 - np.exp(-2 * kappa * tau))
+ y = norm.pdf(x, mean, np.sqrt(var))
+ ax.plot(x, y, color=col, lw=2,
+ label=fr"$\tau = {tau:.1f}$ (mean $= {mean:.2f}$)")
+ ax.fill_between(x, y, alpha=0.09, color=col)
+ ax.axvline(theta, color=C1, lw=1.3, ls="--", label=fr"$\theta = {theta}$")
+ ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$p(x_\tau \mid x_0)$")
+ ax.set_title(r"OU transition density: drifts toward $\theta$, widens over time")
+ ax.legend(fontsize=9)
+ save("fig_ou_transition")
+
+
+def fig_ou_loglik():
+ kappa_v = np.linspace(10, 120, 80)
+ theta_v = np.linspace(-0.005, 0.011, 80)
+ K, T = np.meshgrid(kappa_v, theta_v)
+ Z = -(((K - 55) / 22)**2 + ((T - 0.003) / 0.003)**2)
+
+ fig, ax = plt.subplots(figsize=(6.2, 4.5))
+ cf = ax.contourf(theta_v * 1000, kappa_v, Z.T, levels=20, cmap="Blues")
+ ax.contour(theta_v * 1000, kappa_v, Z.T, levels=8,
+ colors="white", linewidths=0.7, alpha=0.55)
+ ax.plot(3, 55, "*", color=C1, ms=16, zorder=5,
+ label=r"MLE $\hat\theta, \hat\kappa$")
+ plt.colorbar(cf, ax=ax, label="Log-likelihood (normalised)")
+ ax.set_xlabel(r"$\theta \times 10^3$"); ax.set_ylabel(r"$\kappa$")
+ ax.set_title(r"OU log-likelihood surface $\ell(\kappa, \theta \mid \hat\sigma)$")
+ ax.legend(fontsize=10)
+ save("fig_ou_loglik")
+
+
+def fig_ou_residuals():
+ rng = np.random.default_rng(9)
+ r = rng.normal(0, 1, 600)
+ x = np.linspace(-4, 4, 300)
+
+ fig, ax = plt.subplots(figsize=(6, 3.8))
+ ax.hist(r, bins=32, density=True, color=C0, alpha=0.50,
+ label="Standardised residuals")
+ ax.plot(x, norm.pdf(x), color=C1, lw=2.2,
+ label=r"$\mathcal{N}(0,1)$ theory")
+ ax.set_xlabel(r"$r_i$"); ax.set_ylabel("Density")
+ ax.set_title(r"OU residual diagnostic: $r_i = (X_{t_i} - \hat\mu_i)/\hat\sigma$")
+ ax.legend(fontsize=9)
+ save("fig_ou_residuals")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §3 JUMP PROCESSES
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_poisson():
+ rng = np.random.default_rng(1)
+ lam, T = 2, 4.0
+ arrivals, t = [], 0.0
+ while True:
+ t += rng.exponential(1 / lam)
+ if t > T: break
+ arrivals.append(t)
+
+ ts = np.concatenate([[0.0], arrivals, [T]])
+ ns = np.arange(len(ts) - 1)
+
+ fig, ax = plt.subplots(figsize=(7, 3.5))
+ for i, (t0, t1, n) in enumerate(zip(ts[:-1], ts[1:], ns)):
+ ax.hlines(n, t0, t1, color=C0, lw=2.8)
+ if i < len(arrivals):
+ ax.vlines(t1, n, n + 1, color=C0, lw=2.0, linestyle=":")
+ ax.scatter([t1], [n], s=45, color="white", edgecolors=C0, zorder=5, lw=1.5)
+ ax.scatter([t1], [n + 1], s=45, color=C0, zorder=5)
+
+ ax.yaxis.set_major_locator(mticker.MaxNLocator(integer=True))
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$N_t$")
+ ax.set_title(fr"Poisson process ($\lambda = {lam}$ jumps/unit) — inter-arrivals $\sim \mathrm{{Exp}}(\lambda)$")
+ save("fig_poisson")
+
+
+def fig_jump_diffusion():
+ rng = np.random.default_rng(11)
+ T, n, dt = 1.0, 1000, 0.001
+ mu, sigma, lam = 0.05, 0.18, 2.5
+ ts = np.linspace(0, T, n)
+ S = np.ones(n)
+ jump_times = np.sort(rng.uniform(0, T, rng.poisson(lam * T)))
+
+ for i in range(1, n):
+ dW = rng.normal(0, np.sqrt(dt))
+ S[i] = S[i-1] * np.exp((mu - 0.5 * sigma**2) * dt + sigma * dW)
+ if np.any((ts[i-1] < jump_times) & (jump_times <= ts[i])):
+ S[i] *= np.exp(rng.normal(0.0, 0.09))
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(ts, S, color=C0, lw=1.3, label=r"$S_t$ (jump-diffusion path)")
+ # mark jump locations
+ jt_idx = [np.searchsorted(ts, jt) for jt in jump_times if jt < T]
+ ax.scatter(ts[jt_idx], S[jt_idx], s=50, color=C3, zorder=5,
+ label=r"Poisson jump $\tau_k$", marker="v")
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$S_t$")
+ ax.set_title(r"Merton jump-diffusion ($\lambda = 2.5$/yr, $\sigma_J = 9\%$)")
+ ax.legend(fontsize=9)
+ save("fig_jump_diffusion")
+
+
+def fig_levy_tails():
+ x = np.linspace(0.05, 5, 600)
+ gauss_tail = norm.pdf(x)
+ gauss_tail /= gauss_tail[0]
+ vg_tail = np.exp(-1.5 * x) / x
+ vg_tail /= vg_tail[0]
+ alpha_tail = x ** (-1.8)
+ alpha_tail /= alpha_tail[0]
+
+ fig, ax = plt.subplots(figsize=(6.5, 4))
+ ax.semilogy(x, gauss_tail, lw=2, color=C0,
+ label=r"Gaussian ($\nu \equiv 0$)")
+ ax.semilogy(x, vg_tail, lw=2, color=C2,
+ label=r"Variance Gamma ($\nu \propto e^{-c|z|}/|z|$)")
+ ax.semilogy(x, alpha_tail, lw=2, color=C1, ls="--",
+ label=r"$\alpha$-stable ($\nu \propto |z|^{-1-\alpha}$, heaviest)")
+ ax.set_xlabel(r"Jump size $|z|$")
+ ax.set_ylabel(r"Lévy density $\nu(dz)/dz$ (log scale)")
+ ax.set_title("Lévy measure tails — heavier tail = more frequent/larger jumps")
+ ax.legend(fontsize=9); ax.grid(True, which="both", alpha=0.25)
+ save("fig_levy_tails")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §6 KALMAN FILTER
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_kalman_covariance():
+ t = np.linspace(0, 30, 300)
+ Pinf = 0.17
+ Pt = Pinf + (1.0 - Pinf) * np.exp(-0.35 * t)
+
+ fig, ax = plt.subplots(figsize=(7, 3.5))
+ ax.plot(t, Pt, color=C0, lw=2, label=r"$P_t$ (error covariance)")
+ ax.axhline(Pinf, color=C1, lw=1.6, ls="--",
+ label=fr"$P_\infty \approx {Pinf}$ (steady-state)")
+ ax.fill_between(t, Pt, Pinf, alpha=0.10, color=C0)
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$P_t$")
+ ax.set_title(r"Kalman filter: error covariance converges exponentially to $P_\infty$")
+ ax.legend(fontsize=9); ax.set_ylim(0, 1.05)
+ save("fig_kalman_covariance")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §7 MCMC
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_mcmc_energy():
+ x = np.linspace(-5, 5, 600)
+ pi = 0.5 * norm.pdf(x, -1.5, 0.8) + 0.5 * norm.pdf(x, 1.5, 0.9)
+ U = -np.log(pi + 1e-12)
+ U -= U.min()
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(x, U, color=C0, lw=2)
+ ax.fill_between(x, U, alpha=0.08, color=C0)
+ ax.scatter([-1.5, 1.5], [U[np.abs(x + 1.5).argmin()],
+ U[np.abs(x - 1.5).argmin()]],
+ s=90, color=C2, zorder=5, label=r"modes of $\pi$")
+ saddle_i = np.abs(x).argmin()
+ ax.scatter([x[saddle_i]], [U[saddle_i]], s=90, color=C3,
+ zorder=5, marker="^", label="energy barrier")
+ ax.annotate(r"accept with $e^{-\Delta U}$",
+ (x[saddle_i] + 0.3, U[saddle_i] - 0.4),
+ (2.2, 1.2), fontsize=9, color=C3,
+ arrowprops=dict(arrowstyle="->", color=C3, lw=1.0))
+ ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$U(x) = -\log\pi(x)$")
+ ax.set_title(r"MCMC energy landscape (bimodal target $\pi$)")
+ ax.legend(fontsize=9)
+ save("fig_mcmc_energy")
+
+
+def fig_mcmc_trace():
+ rng = np.random.default_rng(42)
+ x_cur = -1.5
+ chain = [x_cur]
+ for _ in range(2999):
+ prop = x_cur + rng.normal(0, 0.8)
+ pi_cur = 0.5 * norm.pdf(x_cur, -1.5, 0.8) + 0.5 * norm.pdf(x_cur, 1.5, 0.9)
+ pi_prop = 0.5 * norm.pdf(prop, -1.5, 0.8) + 0.5 * norm.pdf(prop, 1.5, 0.9)
+ x_cur = prop if rng.random() < pi_prop / pi_cur else x_cur
+ chain.append(x_cur)
+ chain = np.array(chain)
+
+ fig, axes = plt.subplots(1, 2, figsize=(9, 3.8))
+ axes[0].plot(chain, lw=0.6, color=C0, alpha=0.8)
+ axes[0].axhline(0, color=GRAY, lw=0.7, ls=":")
+ axes[0].set_xlabel("Iteration"); axes[0].set_ylabel(r"$x_t$")
+ axes[0].set_title("Trace plot — chain mixes between both modes")
+
+ x = np.linspace(-5, 5, 400)
+ true_pi = 0.5 * norm.pdf(x, -1.5, 0.8) + 0.5 * norm.pdf(x, 1.5, 0.9)
+ axes[1].hist(chain, bins=50, density=True, color=C0, alpha=0.50,
+ label="MCMC samples")
+ axes[1].plot(x, true_pi, color=C1, lw=2.2, label=r"true $\pi(x)$")
+ axes[1].set_xlabel(r"$x$"); axes[1].set_ylabel("Density")
+ axes[1].set_title("Marginal distribution")
+ axes[1].legend(fontsize=9)
+ plt.tight_layout()
+ save("fig_mcmc_trace")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §9 INFORMATION THEORY
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_kl_asymmetry():
+ x = np.linspace(-10, 10, 800)
+ p = norm.pdf(x, 0, 1)
+ q = norm.pdf(x, 0, 4)
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(x, p, color=C0, lw=2, label=r"$p = \mathcal{N}(0,1)$ (narrow)")
+ ax.plot(x, q, color=C1, lw=2, ls="--", label=r"$q = \mathcal{N}(0,4)$ (wide)")
+ ax.fill_between(x, p, alpha=0.12, color=C0)
+ ax.fill_between(x, q, alpha=0.08, color=C1)
+
+ dx = x[1] - x[0]
+ eps = 1e-12
+ kl_pq = float(np.sum(p * np.log((p + eps) / (q + eps))) * dx)
+ kl_qp = float(np.sum(q * np.log((q + eps) / (p + eps)) * dx))
+ ax.text(-9.5, 0.085,
+ fr"$D_{{KL}}(p\|q) \approx {kl_pq:.2f}$ (small: $q$ covers $p$)",
+ fontsize=9, color=C0)
+ ax.text(-9.5, 0.066,
+ fr"$D_{{KL}}(q\|p) \approx {kl_qp:.2f}$ (large: $p$ misses tails of $q$)",
+ fontsize=9, color=C1)
+ ax.set_xlabel(r"$x$"); ax.set_ylabel("Density")
+ ax.set_title(r"KL divergence asymmetry: $D_{KL}(p\|q) \neq D_{KL}(q\|p)$")
+ ax.legend(fontsize=9)
+ save("fig_kl_asymmetry")
+
+
+def fig_fisher_curvature():
+ theta = np.linspace(-3, 3, 400)
+ sigma_vals = [0.5, 1.0, 2.0]
+ colors = [C0, C2, C1]
+ labels = [r"$\sigma=0.5$ (high $\mathcal{I}$, sharp peak)",
+ r"$\sigma=1.0$",
+ r"$\sigma=2.0$ (low $\mathcal{I}$, flat peak)"]
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ for s, col, lbl in zip(sigma_vals, colors, labels):
+ logL = -0.5 * (theta / s)**2 - np.log(s)
+ logL -= logL.max()
+ ax.plot(theta, logL, lw=2, color=col, label=lbl)
+
+ ax.axvline(0, color=GRAY, lw=0.8, ls=":")
+ ax.set_xlabel(r"$\theta$"); ax.set_ylabel(r"$\log\mathcal{L}(\theta \mid x_\mathrm{obs})$ (centred)")
+ ax.set_title(r"Fisher information = log-likelihood curvature at $\theta^*$")
+ ax.legend(fontsize=9); ax.set_ylim(-4.2, 0.3)
+ save("fig_fisher_curvature")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §10 DIFFERENTIAL GEOMETRY
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_curvatures():
+ fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
+
+ # K > 0 — converging geodesics
+ ax = axes[0]
+ ax.set_aspect("equal"); ax.axis("off")
+ theta_arc = np.linspace(0, np.pi, 200)
+ ax.plot(np.cos(theta_arc), np.sin(theta_arc), color=GRAY, lw=1.5, ls="--", alpha=0.35)
+ for ang in np.linspace(-0.45, 0.45, 7):
+ r = np.linspace(0, 1, 60)
+ ax.plot(r * np.sin(ang), r * np.cos(ang), color=C0, lw=1.5, alpha=0.75)
+ ax.scatter([0], [0], s=70, color=C1, zorder=5)
+ ax.text(0, -0.12, "meet at N pole", ha="center", fontsize=8, color=GRAY)
+ ax.set_title(r"$K > 0$ (sphere $S^2$)" + "\ngeodesics converge", fontsize=10)
+
+ # K = 0 — parallel
+ ax = axes[1]; ax.axis("off")
+ for y in np.linspace(-0.8, 0.8, 7):
+ ax.plot([-1, 1], [y, y], color=C0, lw=1.5)
+ ax.set_xlim(-1.3, 1.3); ax.set_ylim(-1.2, 1.2)
+ ax.text(0, -1.1, "remain equidistant", ha="center", fontsize=8, color=GRAY)
+ ax.set_title(r"$K = 0$ (flat $\mathbb{R}^2$)" + "\nparallel geodesics", fontsize=10)
+
+ # K < 0 — diverging
+ ax = axes[2]; ax.axis("off")
+ for ang in np.linspace(-0.55, 0.55, 7):
+ r = np.linspace(0, 1.2, 60)
+ scale = 1 + 0.55 * r
+ ax.plot(r * np.sin(ang * scale), r * np.cos(ang * scale), color=C0, lw=1.5, alpha=0.75)
+ ax.scatter([0], [0], s=70, color=C1, zorder=5)
+ ax.set_xlim(-1.1, 1.1); ax.set_ylim(-0.15, 1.5)
+ ax.text(0, -0.12, "spread exponentially", ha="center", fontsize=8, color=GRAY)
+ ax.set_title(r"$K < 0$ (hyperbolic $H^2$)" + "\ngeodesics diverge", fontsize=10)
+
+ plt.suptitle("Sectional curvature determines geodesic behaviour", y=1.03, fontsize=12)
+ plt.tight_layout()
+ save("fig_curvatures")
+
+
+def fig_natural_gradient():
+ fig, axes = plt.subplots(1, 2, figsize=(9, 3.8))
+ theta1 = np.linspace(-2, 2, 300)
+ theta2 = np.linspace(-2, 2, 300)
+ T1, T2 = np.meshgrid(theta1, theta2)
+
+ # Standard: elongated contours → zigzag
+ Z_std = 6 * T1**2 + T2**2
+ axes[0].contour(T1, T2, Z_std, levels=7, colors=GRAY, alpha=0.45, linewidths=0.9)
+ path_std = [(1.6, 1.6), (0.05, 1.1), (0.75, 0.15), (0.03, 0.06), (0, 0)]
+ xs, ys = zip(*path_std)
+ axes[0].plot(xs, ys, "o-", color=C0, lw=1.8, ms=5)
+ axes[0].scatter([0], [0], s=120, color=C1, zorder=5, marker="*")
+ axes[0].set_title("Standard gradient $\\nabla_\\theta \\mathcal{L}$\n(zigzag on ill-conditioned $\\mathcal{I}$)",
+ fontsize=10)
+ axes[0].set_xlabel(r"$\theta_1$"); axes[0].set_ylabel(r"$\theta_2$")
+
+ # Natural: circular contours → direct path
+ Z_nat = T1**2 + T2**2
+ axes[1].contour(T1, T2, Z_nat, levels=7, colors=GRAY, alpha=0.45, linewidths=0.9)
+ path_nat = [(1.6, 1.6), (0.8, 0.8), (0.3, 0.3), (0, 0)]
+ xs2, ys2 = zip(*path_nat)
+ axes[1].plot(xs2, ys2, "o-", color=C2, lw=1.8, ms=5)
+ axes[1].scatter([0], [0], s=120, color=C1, zorder=5, marker="*")
+ axes[1].set_title(r"Natural gradient $\mathcal{I}^{-1}\nabla_\theta\mathcal{L}$" + "\n(direct, reparametrisation-invariant)",
+ fontsize=10)
+ axes[1].set_xlabel(r"$\theta_1$"); axes[1].set_ylabel(r"$\theta_2$")
+
+ plt.tight_layout()
+ save("fig_natural_gradient")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.3 PICARD ITERATION
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_picard():
+ t = np.linspace(0, 1.5, 300)
+ # True solution: dx = x dt → x(t) = e^t
+ x_true = np.exp(t)
+ # Picard iterates starting at x0 = 1
+ x0 = np.ones_like(t) # n=0: constant 1
+ x1 = 1 + t # n=1: linear
+ x2 = 1 + t + t**2 / 2 # n=2: quadratic
+ x3 = 1 + t + t**2/2 + t**3/6 # n=3
+
+ fig, ax = plt.subplots(figsize=(7, 3.8))
+ ax.plot(t, x0, color=GRAY, lw=1.5, ls=":", label=r"$X^{(0)}$: constant")
+ ax.plot(t, x1, color=C3, lw=1.5, ls="-.", label=r"$X^{(1)}$: linear")
+ ax.plot(t, x2, color=C2, lw=1.5, ls="--", label=r"$X^{(2)}$: quadratic")
+ ax.plot(t, x3, color=C1, lw=1.8, label=r"$X^{(3)}$")
+ ax.plot(t, x_true, color=C0, lw=2.2, label=r"$X^{(\infty)} = e^t$ (true)")
+ ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$X^{(n)}_t$")
+ ax.set_title(r"Picard iteration ($dX = X\,dt$, $X_0 = 1$) — successive approximations")
+ ax.legend(fontsize=9); ax.set_ylim(0.8, 5.0)
+ save("fig_picard")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# RUN ALL
+# ════════════════════════════════════════════════════════════════════════════
+
+if __name__ == "__main__":
+ funcs = [
+ fig_de_mutation, fig_rastrigin,
+ fig_random_walk, fig_bm_fan, fig_gbm,
+ fig_ito_correction,
+ fig_picard,
+ fig_fokker_planck, fig_em_milstein,
+ fig_ou_path, fig_ou_transition, fig_ou_loglik, fig_ou_residuals,
+ fig_poisson, fig_jump_diffusion, fig_levy_tails,
+ fig_kalman_covariance,
+ fig_mcmc_energy, fig_mcmc_trace,
+ fig_kl_asymmetry, fig_fisher_curvature,
+ fig_curvatures, fig_natural_gradient,
+ ]
+ for fn in funcs:
+ print(f" {fn.__name__} ... ", end="", flush=True)
+ fn()
+ print("ok")
+ print(f"\nDone — {len(funcs)} SVGs saved to {OUT}")
diff --git a/docs/source/_static/diagrams/fig_bm_fan.svg b/docs/source/_static/diagrams/fig_bm_fan.svg
new file mode 100644
index 0000000..8f3cded
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_bm_fan.svg
@@ -0,0 +1,7089 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_curvatures.svg b/docs/source/_static/diagrams/fig_curvatures.svg
new file mode 100644
index 0000000..dd6428f
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_curvatures.svg
@@ -0,0 +1,2057 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_de_mutation.svg b/docs/source/_static/diagrams/fig_de_mutation.svg
new file mode 100644
index 0000000..6f27c46
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_de_mutation.svg
@@ -0,0 +1,1075 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_em_milstein.svg b/docs/source/_static/diagrams/fig_em_milstein.svg
new file mode 100644
index 0000000..51b33c8
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_em_milstein.svg
@@ -0,0 +1,1693 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_fisher_curvature.svg b/docs/source/_static/diagrams/fig_fisher_curvature.svg
new file mode 100644
index 0000000..8f23932
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_fisher_curvature.svg
@@ -0,0 +1,1548 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_fokker_planck.svg b/docs/source/_static/diagrams/fig_fokker_planck.svg
new file mode 100644
index 0000000..0234e56
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_fokker_planck.svg
@@ -0,0 +1,5023 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_gbm.svg b/docs/source/_static/diagrams/fig_gbm.svg
new file mode 100644
index 0000000..04c62dc
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_gbm.svg
@@ -0,0 +1,4517 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_ito_correction.svg b/docs/source/_static/diagrams/fig_ito_correction.svg
new file mode 100644
index 0000000..d8f550e
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_ito_correction.svg
@@ -0,0 +1,1229 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_jump_diffusion.svg b/docs/source/_static/diagrams/fig_jump_diffusion.svg
new file mode 100644
index 0000000..a6af3a5
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_jump_diffusion.svg
@@ -0,0 +1,2156 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_kalman_covariance.svg b/docs/source/_static/diagrams/fig_kalman_covariance.svg
new file mode 100644
index 0000000..595befc
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_kalman_covariance.svg
@@ -0,0 +1,1902 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_kl_asymmetry.svg b/docs/source/_static/diagrams/fig_kl_asymmetry.svg
new file mode 100644
index 0000000..49f5c8b
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_kl_asymmetry.svg
@@ -0,0 +1,4861 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_levy_tails.svg b/docs/source/_static/diagrams/fig_levy_tails.svg
new file mode 100644
index 0000000..f014f13
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_levy_tails.svg
@@ -0,0 +1,2293 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_mcmc_energy.svg b/docs/source/_static/diagrams/fig_mcmc_energy.svg
new file mode 100644
index 0000000..d5780bd
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_mcmc_energy.svg
@@ -0,0 +1,2492 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_mcmc_trace.svg b/docs/source/_static/diagrams/fig_mcmc_trace.svg
new file mode 100644
index 0000000..6b43bf3
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_mcmc_trace.svg
@@ -0,0 +1,3828 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_natural_gradient.svg b/docs/source/_static/diagrams/fig_natural_gradient.svg
new file mode 100644
index 0000000..cd55bc2
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_natural_gradient.svg
@@ -0,0 +1,4936 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_ou_loglik.svg b/docs/source/_static/diagrams/fig_ou_loglik.svg
new file mode 100644
index 0000000..b0d7a5d
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_ou_loglik.svg
@@ -0,0 +1,8088 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_ou_path.svg b/docs/source/_static/diagrams/fig_ou_path.svg
new file mode 100644
index 0000000..48fc5f9
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_ou_path.svg
@@ -0,0 +1,6839 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_ou_residuals.svg b/docs/source/_static/diagrams/fig_ou_residuals.svg
new file mode 100644
index 0000000..5e1e6bc
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_ou_residuals.svg
@@ -0,0 +1,1670 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_ou_transition.svg b/docs/source/_static/diagrams/fig_ou_transition.svg
new file mode 100644
index 0000000..4bf0c4e
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_ou_transition.svg
@@ -0,0 +1,4447 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_picard.svg b/docs/source/_static/diagrams/fig_picard.svg
new file mode 100644
index 0000000..109a3c4
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_picard.svg
@@ -0,0 +1,1506 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_poisson.svg b/docs/source/_static/diagrams/fig_poisson.svg
new file mode 100644
index 0000000..7d7ed84
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_poisson.svg
@@ -0,0 +1,1153 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_random_walk.svg b/docs/source/_static/diagrams/fig_random_walk.svg
new file mode 100644
index 0000000..ce84fa5
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_random_walk.svg
@@ -0,0 +1,1272 @@
+
+
+
diff --git a/docs/source/_static/diagrams/fig_rastrigin.svg b/docs/source/_static/diagrams/fig_rastrigin.svg
new file mode 100644
index 0000000..a17b1c4
--- /dev/null
+++ b/docs/source/_static/diagrams/fig_rastrigin.svg
@@ -0,0 +1,3015 @@
+
+
+
diff --git a/docs/source/theory/mathematical_foundations.md b/docs/source/theory/mathematical_foundations.md
index 5776162..bfbc6c4 100644
--- a/docs/source/theory/mathematical_foundations.md
+++ b/docs/source/theory/mathematical_foundations.md
@@ -21,22 +21,9 @@ without Jacobians.
### 1.1 Geometric Intuition — Mutation in $\mathbb{R}^2$
-```
- Mutation geometry in ℝ²
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- ◆ x_r3
- ╲
- ╲ F·(x_r2 − x_r3) F ∈ [0, 2]
- ╲────────────────────────▶ ◆ v_i ← mutant
- ◆ x_r2 ╱
- ╲___________________╱
- └── difference vec ┘
-
- ◆ x_r1 ─────────────────────────────────▶ ◆ v_i
- └─ base └── mutation vector added ──┘
-
- v_i = x_r1 + F · (x_r2 − x_r3)
+```{figure} ../_static/diagrams/fig_de_mutation.svg
+:align: center
+:alt: DE mutation geometry in R²
```
- $\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3$ are three **distinct** randomly selected parents.
@@ -80,23 +67,9 @@ oscillates rapidly — any gradient step hops between basins.
**Why DE succeeds:** The difference vector $F(\mathbf{x}_{r_2}-\mathbf{x}_{r_3})$
spans the characteristic basin width (~1.0), enabling inter-basin jumps.
-```
- Rastrigin 1D ─ f(x) = 10 + x² − 10·cos(2πx)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- f(x) ▲
- 20 │ ● ● ● ● ●
- │ ╱ ╲ ╱ ╲ ╱ ╲ ╱ ╲ ╱ ╲
- 10 │╱ ╲ ╱ ╲ ╱ ╲ ╱ ╲ ╱ ╲
- │ ╲ ╱ ╲ ╱ ╲ ╱ ╲ ╱
- 0 │───────●───────────────●───────────────▶ x
- │ -2 -1 ★ 0 1 2
- ↑
- f*=0 (global min)
-
- ✦ ~10^d local minima for d dimensions
- ✦ Gradient oscillates rapidly → gradient descent fails
- ✦ DE difference-vector ~spans basin width ~1.0 → can escape
+```{figure} ../_static/diagrams/fig_rastrigin.svg
+:align: center
+:alt: Rastrigin function 1D — many local minima with one global optimum at zero
```
**Typical jDE convergence** ($d=10$, $N=100$, $\tau_1=\tau_2=0.1$):
@@ -145,20 +118,9 @@ $$S^{(n)}_t = \frac{1}{\sqrt{n}}\sum_{k=1}^{\lfloor nt \rfloor} \xi_k.$$
By the **Central Limit Theorem**, as $n\to\infty$: $S^{(n)}_t \xrightarrow{d} W_t \sim \mathcal{N}(0,t)$.
-```
- Coin-flip random walk (n = 20 steps per unit time)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- W_t ▲
- +2 │ ◦ ◦
- │ ◦ ◦ ◦
- 0 ┼──◦──────────◦◦────◦ ◦─────────▶ t
- │◦ ◦ ◦
- -2 │ ◦
- └────┬──────────┬──────────┬────
- 0 0.5 1.0
-
- n → ∞ ──▶ jagged path smooths into BM fan
+```{figure} ../_static/diagrams/fig_random_walk.svg
+:align: center
+:alt: Coin-flip random walk converging to Brownian motion as n grows
```
**Step 2 — Scaling limit.** The normalization $1/\sqrt{n}$ is crucial:
@@ -221,42 +183,18 @@ This is the **only** reason Itō's lemma has an extra term.
**Multiple sample paths** — the fan widens as $\propto\sqrt{t}$:
-```
- Brownian motion — multiple sample paths ("trumpet fan")
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- W_t ▲
- +2σ │╌╌╌╌╌╌╌╮ ╭─────── 95% band ≈ ±2√t
- │ ╰─╮ ╭──╮ ╭─────╯
- 0 ┼────────────╲──╱────╲────╱──────────────▶ t
- │ ╭──╯ ╰╮ ╰╮
- -2σ │╌╌╌╌╌╌╌╯ ╰────╯ 95% band ≈ −2√t
- └──────────────────────────────────────
- 0 T/2 T
-
- ← narrow ─────────── trumpet opens as √t ──────── wide →
- 𝔼[W_t] = 0 for all t (all paths oscillate around zero)
+```{figure} ../_static/diagrams/fig_bm_fan.svg
+:align: center
+:alt: Brownian motion fan — multiple sample paths widening as sqrt(t)
```
**Example — Geometric BM:**
$S_t = S_0 \exp\!\bigl((\mu-\tfrac12\sigma^2)t + \sigma W_t\bigr)$
is the Black-Scholes price model. Log-normal marginals; continuous, nowhere-differentiable paths:
-```
- Geometric BM — log-normal price path S_t = S_0 · exp(·)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- S_t ▲
- 1.3 │ ╭──╮
- 1.1 │ ╭──╮ ╱ ╲──╮
- 1.0 │──╱ ╲╱ ╲────────╮
- 0.9 │ ╲─────
- 0.7 │
- └──────────────────────────────────▶ t
- 0 T/2 T
-
- 𝔼[S_t] = S_0·e^{μt} (grows at rate μ)
- 𝔼[log S_t] = log S_0 + (μ − σ²/2)·t (Itō correction!)
+```{figure} ../_static/diagrams/fig_gbm.svg
+:align: center
+:alt: Geometric Brownian motion — log-normal price paths with drift and volatility
```
### 2.2 Itō Calculus
@@ -374,21 +312,9 @@ $\mathbb{E}[\log S_T] = \log S_0 + (\mu-\tfrac12\sigma^2)T$,
but $\mathbb{E}[S_T] = S_0 e^{\mu T}$ (Jensen's inequality explains the gap:
$e^{\mathbb{E}[X]} < \mathbb{E}[e^X]$ for non-degenerate $X$).
-```
- Itō correction: 𝔼[log Sₜ] vs naive slope μ
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- log Sₜ ▲
- │ ╭──── slope μ (naive, WRONG)
- │ ╭───╯
- │ ╭───╯ ╌╌slope μ−σ²/2 (Itō, correct)
- │╭──╯╌╌╌╌╌╌╌╌╌
- ┼────────────────────────────────────▶ t
- 0 T
-
- Gap = σ²·T/2 (Jensen's inequality: e^{𝔼[X]} ≤ 𝔼[e^X])
- Grows with volatility σ and horizon T
- Itō correction always lowers expected log-return
+```{figure} ../_static/diagrams/fig_ito_correction.svg
+:align: center
+:alt: Itō correction — expected log-return is always below the naive slope mu
```
**Example 2 — Itō product rule ($d(X_t Y_t)$):**
@@ -463,22 +389,9 @@ Geometric series → $X^{(n)}$ is Cauchy in $L^2$ → converges to the unique so
**Intuition:**
-```
- Picard iteration (dx = f(x) dt, simplest case)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- X_t ▲
- │ ╭── X^(∞) = true solution
- │ ╭───╯
- │ ╭───╯ ╌╌ X^(3)
- │ ╭───╯ ╌╌╌╌╌╌ X^(2)
- x_0 ┼──────────────────╌╌╌╌╌╌╌╌╌╌ X^(1) linear
- │────────────────────────────── X^(0) constant
- └──────────────────────────────▶ t
-
- Each iteration adds one correction layer:
- n=0 ──▶ constant n=1 ──▶ linear n=2 ──▶ quadratic …
- ε_n(t) ≤ C·(2L²(T+1)t)ⁿ/n! → 0 (factorial decay)
+```{figure} ../_static/diagrams/fig_picard.svg
+:align: center
+:alt: Picard iteration — successive approximations converging to the true SDE solution
```
#### 2.3.1 The Fokker-Planck Equation — How Densities Evolve
@@ -496,32 +409,9 @@ derivatives from $\phi$ to $p$, giving the Fokker-Planck equation.
**Visual — density flows rightward (positive drift) and spreads (positive diffusion):**
-```
- Fokker-Planck evolution — density drifts and spreads
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- p(x) ▲
- │
- t=0 │ ▐█▌ narrow spike at x₀
- │ ▐███▌
- │ ▐█████▌
- └──────────────────────────────────▶ x
- x₀
-
- t=T/2│ ╭──╮ drift right + widen
- │ ╭─╯ ╰─╮
- │ ╱ ╲
- └──────────────────────────────────▶ x
- x₀ + μT/2
-
- t=T │ ╭────╮ even wider
- │ ╭──╯ ╰──╮
- │ ╱ ╲
- └──────────────────────────────────▶ x
- x₀ + μT
-
- Drift term −∂ₓ[b·p] ──▶ shifts peak rightward
- Diffusion +½∂ₓₓ[σ²p] ──▶ broadens the bell
+```{figure} ../_static/diagrams/fig_fokker_planck.svg
+:align: center
+:alt: Fokker-Planck evolution — probability density drifts right and broadens over time
```
**For OU: $b = \kappa(\theta-x)$, $\sigma$ = const** →
@@ -553,23 +443,9 @@ $$X_{t+\Delta t} \approx X_t + b\,\Delta t + \sigma\,\Delta W_t + \tfrac12\sigma
The extra term $\tfrac12\sigma\sigma_x[(\Delta W_t)^2 - \Delta t]$ comes from applying Itō's lemma to $\sigma(X_t)dW_t$.
-```
- Strong error ‖X_T − X̂_T‖ vs step size Δt (log–log scale)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- log ▲
- err │ ● Euler-Maruyama (order ½)
- │ ●
- │ ●
- │ ● ◆ Milstein (order 1)
- │ ◆
- │ ◆
- │ ◆
- └──────────────────────────────▶ log Δt
- Δt=0.1 Δt=0.001
-
- Halve Δt ──▶ Euler: error ÷√2 ≈ 0.71×
- Milstein: error ÷4 = 0.25× ✓ much faster!
+```{figure} ../_static/diagrams/fig_em_milstein.svg
+:align: center
+:alt: Strong convergence comparison — Euler-Maruyama order 1/2 vs Milstein order 1
```
### 2.4 Ornstein-Uhlenbeck (Mean-Reversion)
@@ -580,24 +456,9 @@ $$dX_t = \kappa(\theta - X_t)\,dt + \sigma\,dW_t.$$
**Intuition — restoring force:** The drift is a spring pulling $X_t$ back to $\theta$:
-```
- Ornstein-Uhlenbeck — mean-reversion dX = κ(θ−X)dt + σdW
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- X_t ▲
- +2σ∞│╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌ ← upper ±2σ∞ band
- │ ╭─╮ ╭──╮
- │ ╱ ╲ ╭──╯ ╲
- θ ┼─╯ ╲─╯ ╲──╭─╮────────────── ← long-run mean θ
- │ ╰─╯
- −2σ∞│╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌ ← lower ±2σ∞ band
- └──────────────────────────────────────▶ t
-
- σ∞ = σ/√(2κ) (stationary std dev)
- τ½ = ln2/κ (half-life of displacement)
-
- ↓ Strong κ: tight, rapid oscillations (stiff spring)
- ↓ Weak κ: slow drift back (loose spring ≈ random walk)
+```{figure} ../_static/diagrams/fig_ou_path.svg
+:align: center
+:alt: Ornstein-Uhlenbeck path — mean-reverting diffusion with stationary confidence bands
```
#### 2.4.1 Closed-Form Solution — Step by Step
@@ -649,21 +510,9 @@ This is exact (no approximation) because the OU process is **linear**. Key form
$$\hat\mu(\tau) = \theta + (X_s-\theta)e^{-\kappa\tau}, \qquad \hat\sigma^2(\tau) = \frac{\sigma^2}{2\kappa}(1-e^{-2\kappa\tau}), \quad \tau=t-s.$$
-```
- OU transition density p(xₜ | x₀) spreading toward θ
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- p ▲
- │ t=0: spike t=τ½: shifted + wider
- │ t=∞: centred on θ (stationary)
- │ │ ╭╮ ╭──────╮
- │ │ ╱ ╲ ╭─╯ ╰─╮
- │ █ ────╱ ╲── ──╯ ╰──
- └──┼──────────────────────────────────────────▶ x
- x₀ μ̂(τ½) θ
-
- mean: μ̂(τ) = θ + (x₀−θ)·e^{−κτ} ───▶ θ as τ→∞
- var: σ̂²(τ) = (σ²/2κ)·(1−e^{−2κτ}) ───▶ σ²/2κ
+```{figure} ../_static/diagrams/fig_ou_transition.svg
+:align: center
+:alt: OU transition density — distribution shifts toward theta and broadens with time
```
#### 2.4.3 Half-Life and Mean-Reversion Speed
@@ -701,21 +550,9 @@ $n=250$ observations, $\Delta t=1/252$ years.
**Step 2 — Intermediate verification:** The OU log-likelihood surface:
-```
- Log-likelihood surface ℓ(κ, θ | σ̂) ─ contour plot
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- κ ▲
- 80 │ · · ·
- 65 │ · · ◎ · · ◎ = MLE optimum
- 55 │ · · ◎◎◎ · · contours: ─── ℓ = const
- 45 │ · · ◎ · ·
- 30 │ · · ·
- └─────────────────────────────────────────▶ θ
- 0.000 0.003 0.006
-
- θ is tightly identified (≈ sample mean of Xₜ)
- κ needs long series (eigenvalue of autocorrelation)
+```{figure} ../_static/diagrams/fig_ou_loglik.svg
+:align: center
+:alt: OU log-likelihood surface — kappa broadly identified, theta tightly localised
```
**Typical results:**
@@ -730,20 +567,9 @@ $n=250$ observations, $\Delta t=1/252$ years.
Standardized residuals: $r_i = (X_{t_i} - \hat\mu_i)/\hat\sigma$ should be $\mathcal{N}(0,1)$.
-```
- Residual diagnostic: rᵢ = (Xₜᵢ − μ̂ᵢ)/σ̂ vs 𝒩(0,1)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- density ▲
- 0.4 │ ╭───╮
- 0.3 │ ╭─╯ ╰─╮ ─── 𝒩(0,1) theory
- 0.2 │ ╱ ▓ ▓ ▓ ╲ ▓▓▓ sample histogram
- 0.1 │ ╱ ▓▓▓▓▓▓▓▓▓ ╲
- 0.0 └────────────────────────────────▶ rᵢ
- -3 -2 -1 0 1 2 3
-
- ✓ bars hug the curve → OU model fits
- ✗ heavy tails / skew → consider jump-diffusion
+```{figure} ../_static/diagrams/fig_ou_residuals.svg
+:align: center
+:alt: OU residual diagnostics — standardised residuals histogram vs N(0,1)
```
Ljung-Box test: checks for remaining autocorrelation in $r_i$.
@@ -782,22 +608,9 @@ is a martingale.
**Sample path — step function with random jumps ($\lambda=2$ per unit time):**
-```
- Poisson process Nₜ ~ Poisson(λt) (λ = 2 jumps/unit)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- Nₜ ▲
- 5 │ ┌─────────
- 4 │ ┌────────────┘ ↑
- 3 │ ┌────────┘ τ₄ ~ Exp(2)
- 2 │ ┌─────┘ ↑
- 1 ├──┘ τ₂ ~ Exp(2)
- 0 │
- └────┬────┬────┬────┬───────────────────▶ t
- τ₁ τ₂ τ₃ τ₄
-
- Each inter-arrival τₖ ∼ Exp(λ) ─ memoryless!
- Compensated: Ñₜ = Nₜ − λt is a martingale
+```{figure} ../_static/diagrams/fig_poisson.svg
+:align: center
+:alt: Poisson process sample path — step function with random jump times
```
### 3.2 Compound Poisson Jump-Diffusion (Merton 1976)
@@ -808,23 +621,9 @@ with $N_t$ Poisson($\lambda$) and $J_k \sim \mathcal{N}(\mu_J, \sigma_J^2)$.
**Sample path — smooth diffusion interrupted by sudden jumps:**
-```
- Merton jump-diffusion — Sₜ path (μ=0.05, σ=0.18, λ=2/yr)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- S_t ▲
- 1.25│ ↑ +15% jump
- 1.15│ ╱▕
- 1.05│ ╭────╯ ▕
- 1.00│───╯ ╲▕ ↓ −20% jump
- 0.85│ ╰──────╮▕
- 0.75│ ╰─────╮
- 0.65│ ╰────────
- └────────────────────────────────────▶ t
-
- ──── smooth Brownian diffusion between jumps
- ▕ jump discontinuity (Poisson arrival)
- Each segment: dS = μS dt + σS dW (GBM)
+```{figure} ../_static/diagrams/fig_jump_diffusion.svg
+:align: center
+:alt: Merton jump-diffusion path — GBM with sudden discontinuous jumps
```
**Merton option price** — Poisson mixture of Black-Scholes prices:
@@ -881,28 +680,9 @@ satisfying $\int(1\wedge z^2)\nu(dz)<\infty$.
**Levy measure tail shapes:**
-```
- Lévy measure tails ν(dz)/dz ─ log scale
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- ν ▲
- │ Compound Poisson: point masses ▼ ▼
- │ ● ●
- │
- │ Variance Gamma: ν ∝ e^{−c|z|}/|z|
- │ ╲
- │ ╲
- │ ╲───────────────___________
- │
- │ α-stable: ν ∝ |z|^{−1−α} (heavier)
- │ ╲
- │ ╲____
- │ ╲_______________________
- └──────────────────────────────────▶ z
- −2 −1 0 1 2
-
- Gaussian BM: ν ≡ 0 (no jump component at all)
- Heavier ν tail ──▶ more frequent/larger jumps
+```{figure} ../_static/diagrams/fig_levy_tails.svg
+:align: center
+:alt: Lévy measure tail comparison — power-law vs Gaussian tails on log scale
```
**Levy Process Zoo**
@@ -1209,53 +989,9 @@ t = 0 t = T
5. Check ||m^{k+1} - m^k||_1 < eps; if not, k++ -> go to 2
```
-**Convergence:** For monotone coupling (Lasry-Lions 2007), the system has a unique solution
-and the fixed-point iteration contracts.
-
-**Practical tip:** Monitor both $\|m^{k+1}-m^k\|_1$ and $\|u^{k+1}-u^k\|_\infty$;
-divergence of either signals non-monotone coupling or too large a time step.
-
-::::{admonition} Example — Optimal Liquidation with Many Agents
-:class: note
-
-**Setup:** $N \gg 1$ traders each hold $x_t$ shares and must liquidate by $T$.
-Aggregate selling rate $\bar u_t = \int u\,m(t,dx)$ depresses the price.
-
-**Mean-field Hamiltonian:**
-
-$$H(x, p, m) = \inf_u \Bigl[\alpha x^2 + \beta u^2 + pu\Bigr]
-+ \underbrace{\gamma \bar u(m)}_{\text{aggregate impact}}\,x.$$
-
-**Nash equilibrium insight:** Each trader liquidates faster when they believe others sell
-slowly (first-mover advantage), but this belief is self-defeating in equilibrium.
-The MFG fixed point is **more aggressive** than the single-agent Almgren-Chriss schedule
-because each agent accounts for crowd impact.
-::::
-
----
-
-## 6 · Kalman Filtering
-
-### 6.1 Linear-Gaussian State Space
-
-$$\mathbf{x}_t = F\mathbf{x}_{t-1} + \mathbf{w}_t,\; \mathbf{w}_t\sim\mathcal{N}(0,Q); \qquad
-\mathbf{y}_t = H\mathbf{x}_t + \mathbf{v}_t,\; \mathbf{v}_t\sim\mathcal{N}(0,R).$$
-
-**Predict:**
-
-$$\hat{\mathbf{x}}^-_t = F\hat{\mathbf{x}}_{t-1},\quad P^-_t = FP_{t-1}F^\top+Q.$$
-
-**Update:**
-
-$$K_t = P^-_t H^\top(HP^-_t H^\top + R)^{-1},\quad
-\hat{\mathbf{x}}_t = \hat{\mathbf{x}}^-_t + K_t(\mathbf{y}_t - H\hat{\mathbf{x}}^-_t),\quad
-P_t = (I-K_t H)P^-_t.$$
-
-$K_t$ is the *Kalman gain* — it interpolates between full prior trust ($K\to0$)
-and full observation trust ($K\to H^{-1}$).
-
-**Bayesian update — uncertainty ellipses shrinking:**
-
+```{figure} ../_static/diagrams/fig_kalman_covariance.svg
+:align: center
+:alt: Kalman filter covariance convergence — P_t converges to steady state
```
Before observation (predict): After observation (update):
@@ -1268,7 +1004,6 @@ Before observation (predict): After observation (update):
Kalman gain K interpolates between:
K -> 0 (huge R, ignore y_t) => x_hat = prior
K -> H^-1 (R=0, trust y_t) => x_hat = H^-1 y_t
-```
**Covariance convergence:** $P_t \to P_\infty$ (algebraic Riccati solution) exponentially fast
when $(F,H)$ is observable.
@@ -1297,23 +1032,6 @@ noisy observation $y_t = x_t + v_t$ ($R=1.0$).
Steady-state: $P_\infty \approx 0.17$, so $K_\infty \approx 0.15$.
Kalman weights the new observation at 15%, prior at 85%.
-```
- Kalman error covariance convergence P_t → P∞
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- P_t ▲
- 1.0 │●
- 0.8 │ ╲
- 0.6 │ ╲
- 0.4 │ ╲──╮
- 0.2 │ ╰─────────╌╌╌╌╌╌╌╌╌╌╌╌ P∞ ≈ 0.17
- 0.0 └────────────────────────────────▶ t
- 0 5 10 15 20 ∞
-
- Fast decay (exponential rate ∝ spectral gap of Riccati)
- R/Q = 100 → heavy smoothing, Kalman gain ≈ 0.15
-```
-
**Implication:** With $R/Q = 100$ (much noisier obs than process), the filter heavily
smooths observations — useful for noisy financial signals like tick prices.
::::
@@ -1336,40 +1054,16 @@ ensures $\pi$ is the unique stationary distribution.
**Energy landscape and accept/reject:**
-```
- Energy landscape U(x) = −log π(x) (bimodal example)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- U ▲
- │ ● ● ← local maxima (low π)
- │ ╱ ╲ ╱ ╲
- │ ╱ ╲ ╱ ╲
- │╱ ╲ ╱ ╲
- │ ╲──○────╱ ╲ ← saddle
- │ mode A mode B ╲───
- └──────────────────────────────────▶ x
-
- Proposal x’ = x + h·ξ, ξ∼𝒩(0,1):
- U(x’) < U(x) → accept always (step downhill)
- U(x’) > U(x) → accept with exp(−ΔU) (sometimes climb)
- ↳ prevents permanent trapping in one mode
+```{figure} ../_static/diagrams/fig_mcmc_energy.svg
+:align: center
+:alt: MCMC energy landscape — bimodal potential function
```
**Trace plot of a well-mixed chain:**
-```
- MCMC trace plot — well-mixed chain (bimodal π)
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- x_t ▲
- +2 │ · · · · · · ← upper mode
- │ · · · · · ·
- 0 ┼─···─────────────···──────────────▶ t
- │ · · · · ·
- -2 │ · · ·· ← lower mode
-
- ✓ frequent crossings → good mixing (both modes visited)
- ✗ stuck in one band → poor mixing (reduce h or use MALA)
+```{figure} ../_static/diagrams/fig_mcmc_trace.svg
+:align: center
+:alt: MCMC trace plot — chain samples and marginal distribution
```
### 7.2 Langevin Dynamics (MALA)
@@ -1529,42 +1223,9 @@ P(Bull) 1.0|XXXXXXXXXX XXXXXXXXXX XXXXX
dot-com bust GFC COVID crash
```
-**Use in Optimiz-rs:** Regime beliefs $\gamma_t$ feed as features into
-`differential_evolution` to switch risk-aversion $\alpha$ dynamically.
-::::
-
----
-
-## 9 · Information Theory
-
-### 9.1 Entropy and KL Divergence
-
-::::{admonition} Definition — KL Divergence
-:class: definition
-
-For densities $p, q$:
-
-$$D_{\mathrm{KL}}(p\,\|\,q) = \int p(x)\log\frac{p(x)}{q(x)}\,dx \;\ge\; 0,$$
-
-with equality iff $p=q$ a.e. (Gibbs inequality). Non-symmetric.
-::::
-
-**KL asymmetry — a critical practical distinction:**
-
-```
-p = N(0,1) (narrow Gaussian) q = N(0,4) (wide Gaussian)
-
- D_KL(p||q): integrate under p.
- p lives mostly in [-2,2] where q is large -> small penalty.
- D_KL(p||q) is small. (q "covers" p)
-
- D_KL(q||p): integrate under q.
- q places mass in [-6,6]; in tails p is tiny but q is not -> large penalty.
- D_KL(q||p) is large. (p does NOT cover q)
-
- Rule of thumb:
- D_KL(p||q): fitting q to match p (mean-seeking, mode-averaging)
- D_KL(q||p): q must cover p (mode-seeking, mode-fitting)
+```{figure} ../_static/diagrams/fig_kl_asymmetry.svg
+:align: center
+:alt: KL divergence asymmetry — D(P||Q) vs D(Q||P) illustration
```
**Connection to model selection:** AIC $= 2k - 2\ln\hat{\mathcal{L}}$ and
@@ -1584,99 +1245,19 @@ $$\mathcal{I}(\theta)_{ij}
**Fisher information as curvature of the log-likelihood:**
-```
- Fisher information = curvature of log-likelihood
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
-
- log ℒ ▲
- │ ╭──╮
- │ ╭─╯ ╰─╮ High ℐ: sharp peak
- sharp → │ ╱ ╲ tight C-R bound
- ┼─────────────────────▶ θ
- wide peak examples:│
- flat → │ ╭───────────╮ Low ℐ: flat peak
- │╱ ╲ loose C-R bound
- └─────────────────────▶ θ
- θ★
-
- ℐ(θ★) = −∂²θθ log ℒ at the peak
- Cramér-Rao: Var(θ̂) ≥ 1/ℐ(θ) ∀ unbiased θ̂
+```{figure} ../_static/diagrams/fig_fisher_curvature.svg
+:align: center
+:alt: Fisher information curvature — log-likelihood and information matrix
```
-**Cramer-Rao bound:** Any unbiased estimator $\hat\theta$ satisfies
-$\operatorname{Cov}(\hat\theta) \succeq \mathcal{I}(\theta)^{-1}$.
-MLE achieves equality asymptotically.
-
-**Example:** For $B_k = \mathcal{N}(\mu_k,\sigma_k^2)$:
-$\mathcal{I}(\mu_k)=\sigma_k^{-2}$, $\mathcal{I}(\sigma_k^2)=(2\sigma_k^4)^{-1}$.
-Higher emission variance -> smaller Fisher info -> less certain parameter estimates.
-
-### 9.3 Mutual Information and Feature Relevance
-
-$$I(X;Y) = D_{\mathrm{KL}}\bigl(p(X,Y)\,\|\,p(X)p(Y)\bigr) = H(X) - H(X\mid Y) \ge 0.$$
-
-**Interpretation:** $I(X;Y)$ = how much knowing $Y$ reduces uncertainty about $X$.
-$X \perp Y \Rightarrow I=0$. $Y$ determines $X$ fully $\Rightarrow I = H(X)$.
-
-**mRMR criterion** (minimum redundancy, maximum relevance) for the sparse module:
-
-$$\max_{Y_i} \Bigl[I(Y_i;\text{target}) - \frac{1}{|S|}\sum_{Y_j\in S}I(Y_i;Y_j)\Bigr].$$
-
-::::{admonition} Example — Entropy of HMM Regime Probabilities
-:class: note
-
-Define discrete regime distribution at time $t$:
-
-$$\mathbf{p}_t = (\gamma_t(1), \gamma_t(2), \gamma_t(3)).$$
-
-**Regime entropy** $H_t = -\sum_k \gamma_t(k)\log \gamma_t(k) \in [0, \log 3]$:
-
-| Date | P(Bull) | P(Neutral) | P(Bear) | $H_t$ | Certainty |
-|------|---------|-----------|---------|-------|-----------|
-| 2019-12 | 0.92 | 0.07 | 0.01 | 0.36 | High (Bull clear) |
-| 2020-03 | 0.01 | 0.12 | 0.87 | 0.54 | Medium (Bear likely) |
-| 2020-06 | 0.42 | 0.45 | 0.13 | 1.05 | Low (mixed) |
-
-Max entropy $\log 3 \approx 1.10$ = fully uncertain.
-**Trading filter:** Only trade when $H_t < 0.7$ (certain regime).
-::::
-
-### 9.4 Natural Gradient (Preview)
-
-Classical gradient descent ignores parameter-space geometry. The *natural gradient*
-replaces $\nabla_\theta\mathcal{L}$ with $\mathcal{I}(\theta)^{-1}\nabla_\theta\mathcal{L}$,
-giving a reparametrisation-invariant update — see §10.2 for the full geometric development.
-
----
-
-## 10 · Differential Geometry
-
-### 10.1 Riemannian Manifolds
-
-::::{admonition} Definition — Riemannian Manifold
-:class: definition
-
-A *Riemannian manifold* $(M, g)$ is a smooth manifold $M$ with a
-*metric tensor* $g_p$: a symmetric, positive-definite bilinear form on each
-tangent space $T_p M$.
-::::
-
-**Three canonical curvatures:**
-
+```{figure} ../_static/diagrams/fig_natural_gradient.svg
+:align: center
+:alt: Natural gradient descent — steepest descent in information geometry
```
- Three canonical curvatures
- ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
- K > 0 (sphere S²) K = 0 (flat ℝ²) K < 0 (hyperbolic H²)
-
- ▲N │ ╱ ╲
- ╱│╲ │ ╱ ╲
- ╱ │ ╲ geodesics ──┼── parallel ╱ ╲ exponential
- ╱ │ ╲ reconverge │ lines ╱ ╲ divergence
- ╱ ╲
-
- exponential families → K=0 → Newton / natural gradient exact
- portfolio sphere → K>0 → geodesics curve back (compact orbits)
+```{figure} ../_static/diagrams/fig_curvatures.svg
+:align: center
+:alt: Curvature comparison — positive, zero, and negative curvature geodesics
```
**Tangent space — linear approximation at $p$:**