ML Is Bayesian, and Gen ML Requires Bayes

Code

theme = ({
  blue: "#2f64ad", amber: "#c2680c", slate: "#64748b", ink: "#334155",
  blueFill: "rgba(47,100,173,0.14)", amberFill: "rgba(194,104,12,0.12)",
  blues4: ["#b9d0ea", "#7ea9d8", "#4c80c2", "#2f64ad"]
})
themedPlot = (o = {}) => Plot.plot({
  ...o,
  style: {fontFamily: "inherit", fontSize: "12.5px", background: "transparent",
          overflow: "visible", ...(o.style || {})},
  x: {tickSize: 0, ...(o.x || {})},
  y: {tickSize: 0, ...(o.y || {})}
})

By the end of this page you will be able to look at the entire modern generative toolbox - VAEs, diffusion models, GPT - and see a catalogue of responses to one integral. Better: you’ll have derived the shape of each response yourself, before meeting a single neural network, using nothing fancier than ridge regression and the chain rule of probability. The supervised material that fills the first two thirds of this page is not a warm-up. It is a spring being wound. Watch for three questions that supervised learning gets to dodge. Generative modeling doesn’t dodge anything.

Prerequisites. Comfort with probability at the level of a first course. If you know what a conditional distribution is and you believe the chain rule, you have almost everything you need - which is, quietly, this page’s entire thesis.

1. A model is a distribution

Here is a reframe that costs nothing and changes everything. A regression model is usually introduced as a function: feed in \(x\), get out a best guess \(\hat y = f(x)\). But look at what you actually did when you fit it. You minimized squared error - and minimizing squared error is exactly maximizing the likelihood of a Gaussian: \(p_\theta(y \mid x) = \mathcal N(y \mid f_\theta(x), \sigma^2)\). Classification with cross-entropy? A categorical likelihood. You were never fitting functions. You were fitting conditional distributions with the distribution hidden inside a loss.

The standard treatments: MacKay, Information Theory, Inference, and Learning Algorithms (free online); Bishop, Pattern Recognition and Machine Learning.

Why care? Because a function can answer one question - what’s your best guess? - while a distribution answers that one plus three more: how sure are you?, what else is plausible?, and, the sleeper, give me a new sample. Hold on to that last one. Nothing in this section needs it yet. That’s precisely the point.

Figure 1: The same fitted model, described two ways. Left: a function - one guess per x. Right: what it actually is - a conditional density at every x (three of them drawn), with the band as the region holding 95% of each.

2. The fundamental problem

Two distributions sit at the center of everything.

On one side, the empirical distribution \(\hat p\): delta spikes of mass \(1/N\) sitting on your training points. It is the truth as far as you know - every observation you have, and nothing you don’t. On the other side, the model distribution \(p_\theta\): your smooth, parameterized approximation. Learning is the act of pushing the second toward the first.

Which forces the core question of the entire field: how do we measure how close two distributions are?

The honest answer is the Kullback–Leibler divergence:

\[ \mathrm{KL}(P \,\|\, Q) \;=\; \mathbb E_{x \sim P}\!\left[\log \frac{P(x)}{Q(x)}\right]. \]

Read it as an expected surprise penalty: if the world runs on \(P\) and you insist on describing it with \(Q\), KL is how much extra surprise you eat, on average, for your stubbornness. It is not symmetric - \(\mathrm{KL}(P\|Q) \neq \mathrm{KL}(Q\|P)\) - and the asymmetry is not a defect. The two directions penalize different sins, and you should feel the difference in your hands before you trust it.

Code

viewof klDir = Inputs.radio(["KL(P‖Q): mass-covering", "KL(Q‖P): mode-seeking"], {value: "KL(P‖Q): mass-covering", label: "Direction"})
viewof klMu = Inputs.range([-4, 4], {value: 2.4, step: 0.05, label: "Q mean"})
viewof klSigma = Inputs.range([0.3, 3], {value: 0.7, step: 0.05, label: "Q width"})

Figure 2: P is the fixed two-humped density; Q is yours. Drag Q around and switch directions. KL(P‖Q) explodes wherever P has mass that Q ignores - so minimizing it stretches Q to cover both humps. KL(Q‖P) only cares where Q puts mass - so minimizing it lets Q pick one hump and commit. You will meet the symmetrized cousin of this widget (Jensen–Shannon) when we get to GANs.

Figure 3: The same experiment, run exhaustively so the behavior is on the record: each KL direction as a function of Q’s mean, at four widths. Left: KL(P‖Q) is punished brutally whenever Q’s mass fails to cover a hump - note how narrow Qs (small σ) pay enormous penalties everywhere except dead-center on a mode, and how the only way to be safe everywhere is to be wide. Right: KL(Q‖P) is calm by comparison - a narrow Q parked on either hump scores near zero, and the penalty spikes only where Q itself puts mass in P’s valleys. Covering versus committing, as fixed curves.

Now the result that turns KL from philosophy into practice. Take \(P = \hat p\) (the empirical spikes) and \(Q = p_\theta\), expand, and watch the terms fall:

The algebra, for the willing

\[ \mathrm{KL}(\hat p \,\|\, p_\theta) = \underbrace{\sum_i \tfrac1N \log \tfrac{1/N}{1}}_{\text{constant in } \theta} \;-\; \frac1N \sum_{i=1}^N \log p_\theta(x_i). \] The first term is the (fixed) entropy of the empirical distribution. Minimizing KL over \(\theta\) therefore ignores it entirely, leaving only the negative mean log-likelihood.

Minimizing \(\mathrm{KL}(\hat p \| p_\theta)\) is maximizing likelihood. MLE was never a recipe someone handed down; it is distance-minimization toward the empirical distribution, with KL as the ruler. Every loss function you have ever minimized was this, wearing a costume.

3. Under pressure: the collapse

And now the trap springs.

If the target of maximum likelihood is the empirical distribution, then ask: what is the best possible model, the global optimum of the objective, when nothing constrains us? It is the empirical distribution itself. Delta spikes on the training data. Infinite density at every point you observed, zero everywhere else, a training likelihood of \(+\infty\) - and absolutely nothing to say about any point you haven’t seen.

This is worth sitting with, because the standard vocabulary undersells it. We say flexible models “can overfit,” as if overfitting were an occasional hazard. No: unconstrained maximum likelihood aims at the spikes. Overfitting is not the failure mode. It is the optimum.

Do it yourself.

Code

viewof cmBw = Inputs.range([-2, 0.4], {value: -0.2, step: 0.02, label: "log₁₀ bandwidth σ", transform: x => x})
viewof cmTest = Inputs.toggle({label: "Reveal held-out points", value: false})

Figure 4: Twelve training points and a kernel density whose bandwidth you control. Drag σ down and watch the training log-likelihood scream toward +∞ as the density collapses onto spikes - then reveal the held-out points and watch their likelihood crater. You are performing mass collapse: all probability on what we saw, none anywhere else.

Give the pathology its name - mass collapse - because you will meet it again wearing different clothes. (When we get to autoencoders, a deterministic encoder commits exactly this crime in latent space: a Dirac spike of posterior mass per training point, nothing in between. Same disease, new organ.)

Here is the first spring being wound. Supervised learning gets to paper over collapse: hold out a validation set, stop training early, and - the quiet cheat - only ever query the model near data it has seen. A generative model has no such luxuries. Samples come from wherever the mass is. If the mass is wrong, the output is garbage, and no validation-set etiquette will save you. File that away.

4. Beliefs or degeneracy

So constraint is not optional. But before reaching for a fix, be honest about what the collapse means, because the honesty is where Bayes comes from.

Collapse is not irrational. Collapse is rational under zero information. If you genuinely believe nothing about the world beyond the twelve points you observed - no smoothness, no boundedness, no “tomorrow resembles today” - then the spikes are the correct conclusion. The absurdity of that conclusion is a proof by contradiction: you do believe something. You always did. The only question on the table is whether you will say your beliefs out loud, where they can be examined, or leave them implicit, where they can’t.

Bayes’ rule is the bookkeeping for saying them out loud. Set the notation once, cleanly, because it carries the rest of the course. Let \(\theta\) be the set of unknowns we want to learn about - coefficients, weights, latent structure, whatever the model hides from us. Let \(\mathcal D\) be the set of knowns: the data, observed \(N\) times. Then:

\[ \underbrace{P(\theta \mid \mathcal D)}_{\text{posterior}} \;=\; \frac{\overbrace{P(\mathcal D \mid \theta)}^{\text{likelihood}}\;\;\overbrace{P(\theta)}^{\text{prior}}}{\underbrace{P(\mathcal D)}_{\text{normalizing constant}}}. \]

Three named parts. The likelihood is old news - it is the objective we’ve been maximizing since Section 2, now recognized for what it is. The prior is the new arrival: your beliefs about \(\theta\), stated as a distribution, before the data get a word in. And the normalizing constant \(P(\mathcal D)\) is the part we get to ignore for now - it doesn’t involve \(\theta\), so it cannot move any maximization over \(\theta\), and we will happily write \(\propto\) and pretend it isn’t there. Ignore it. But do not forget it. It is lurking in the corner of this entire page like Ghostface, and when it finally strikes - Section 7 - it will turn out to be the most important number in the course.

Take logs and negate, and the posterior becomes an energy with two terms:

\[ -\log P(\theta \mid \mathcal D) \;=\; \underbrace{-\log P(\mathcal D \mid \theta)}_{\text{loss}} \;+\; \underbrace{-\log P(\theta)}_{\text{regularizer}} \;+\; \text{const.} \]

Minimizing that energy has a name. The maximum a posteriori estimate is the peak of the posterior:

\[ \hat\theta_{\mathrm{MAP}} \;=\; \arg\max_\theta \; P(\theta \mid \mathcal D) \;=\; \arg\min_\theta \; \Big[ -\log P(\mathcal D \mid \theta) \;-\; \log P(\theta) \Big], \]

where the normalizing constant drops out of the argmax entirely - Ghostface naps on. MAP is why Bayesian thinking sneaks into practice so easily: it converts inference into optimization - one point, one objective, gradient descent’s home turf - and it is literally the computation every regularized fit performs, whether its author would use the word or not. But hold the concept loosely, because we are moving quickly toward territory where it stops being enough. An argmax is a point, and a point is not a distribution: it has a location and no mass, a peak and no shape. When we reach generative modeling at the end of this page, that deficiency stops being philosophical - generation means drawing from a distribution, and no amount of cleverness makes a single point drawable. Generative ML cannot run on MAP alone. Remember that sentence; the final section is built on it.

Loss plus regularizer. Every regularized objective you have ever minimized was a log-posterior with the prior left unnamed. This identity - call it the Bayesian Reveal - is the keystone of everything downstream: not that you should use priors, but that you already are, and have been since your first homework. And you shouldn’t take a claim like that on vibes. Here are the two canonical unmaskings, worked in full:

Proof: the ℓ₂ penalty is a Gaussian prior (Ridge, unmasked)

Model: \(y_i = \alpha + \mathbf x_i^\top \boldsymbol\beta + \varepsilon_i\) with \(\varepsilon_i \sim \mathcal N(0,\sigma^2)\), independent priors \(\beta_j \sim \mathcal N(0, \tau^2)\) on the slopes, and no prior on the intercept \(\alpha\). (Leaving \(\alpha\) unshrunk is the standard convention, and honestly it’s more pathology than principle: the intercept mostly recenters predictions to the scale of your measurement units, and “the intercept is probably near zero” is a belief about units, not about the world. So we decline to state a belief we don’t hold.)

The negative log posterior, dropping \(\theta\)-free constants:

\[ -\log P(\alpha,\boldsymbol\beta \mid \mathcal D) = \frac{1}{2\sigma^2}\sum_{i=1}^N \big(y_i - \alpha - \mathbf x_i^\top\boldsymbol\beta\big)^2 + \frac{1}{2\tau^2}\sum_{j=1}^P \beta_j^2 + \text{const}. \]

Multiply through by \(2\sigma^2\) - a positive constant, so the argmin doesn’t move:

\[ \sum_{i=1}^N \big(y_i - \alpha - \mathbf x_i^\top\boldsymbol\beta\big)^2 \;+\; \lambda \sum_{j=1}^P \beta_j^2, \qquad \lambda = \frac{\sigma^2}{\tau^2}. \]

Ridge regression. Exactly. The tuning parameter you’ve been cross-validating is a ratio of variances: how noisy you believe the world is, over how large you believe the coefficients are.

Proof: the ℓ₁ penalty is a Laplace prior (Lasso, unmasked)

Same model, same unshrunk intercept; now put \(\beta_j \sim \mathrm{Laplace}(0, b)\) on the slopes, with density \(P(\beta_j) = \frac{1}{2b}\exp\!\big(-|\beta_j|/b\big)\):

\[ -\log P(\alpha,\boldsymbol\beta \mid \mathcal D) = \frac{1}{2\sigma^2}\sum_{i=1}^N \big(y_i - \alpha - \mathbf x_i^\top\boldsymbol\beta\big)^2 + \frac{1}{b}\sum_{j=1}^P |\beta_j| + \text{const}. \]

Multiply through by \(2\sigma^2\):

\[ \sum_{i=1}^N \big(y_i - \alpha - \mathbf x_i^\top\boldsymbol\beta\big)^2 \;+\; \lambda \sum_{j=1}^P |\beta_j|, \qquad \lambda = \frac{2\sigma^2}{b}. \]

The lasso. Exactly. And the famous sparsity of \(\ell_1\) solutions stops being folklore: it’s the Laplace density’s sharp peak at zero doing what sharp peaks do to posteriors.

Jaynes, Probability Theory: The Logic of Science, for the philosophical frame; Murphy, Probabilistic Machine Learning, for the identities. The reveal is mathematically uncontroversial - what’s uncommon is taking it seriously as an organizing lens.

From priors on parameters to priors on functions. Now the observation that takes the blindfold off. In practice, machine learning essentially never moves the prior mean off zero. Nobody centers their weight decay at \(\beta = 0.7\); nobody has ever cross-validated the prior’s location. So if the mean is always zero, ask the uncomfortable question: what is the choice of variance - the choice of penalty family - actually doing? It cannot be encoding where the coefficients are. It must be encoding something else.

It is. Push the prior on parameters through the model, and it becomes a prior on functions: a probability distribution over the explanations you consider believable before any data arrive. Make it concrete. Suppose the truth is quadratic, but you hedge your bets and fit a degree-9 polynomial - capacity to spare, as modern ML always carries. The prior on those nine coefficients is secretly a bet about which of the degree-9 family’s inhabitants are plausible. Draw the whole distribution and look:

Code

// shared data for the fan + the TV density (computed once; only the overlay resample button is reactive)
fanData = {
  // Variance fixed at 1: f is linear in β, so scaling the matched variance
  // rescales every quantile and every TV by √v identically - a slider here
  // is pure zoom, not information. (Confirmed by trying it.)
  const v = 1;
  const tau = Math.sqrt(v), b = Math.sqrt(v / 2);
  const deg = 9, S = 8000;
  const xs = linspace(-1, 1, 61);
  // dark-blue core → light-blue fringe [restyle hexes to site palette later]
  const levels = [
    {p: 0.99, c: "#dbeafe"}, {p: 0.90, c: "#a9c9ef"},
    {p: 0.75, c: "#6b9fdc"}, {p: 0.50, c: "#2f64ad"}
  ];
  const lapDraw = (rng, bb) => {
    const u = rng() - 0.5;
    return -bb * Math.sign(u) * Math.log(Math.max(1 - 2 * Math.abs(u), 1e-12));
  };
  const evalPoly = (beta, x) => {
    let y = 0, xp = 1;
    for (let k = 0; k < deg; k++) { xp *= x; y += beta[k] * xp; }
    return y;
  };
  const out = {v, tau, b, deg, xs, levels, bands: [], tv: [], yext: 0, lapDraw, evalPoly};
  for (const prior of ["Gaussian prior (ℓ₂)", "Laplace prior (ℓ₁)"]) {
    const g = prior.startsWith("Gauss");
    const rngB = mulberry32(g ? 1111 : 2222);       // fixed seeds: stable fan
    const cols = xs.map(() => []);
    for (let sI = 0; sI < S; sI++) {
      const bt = Array.from({length: deg}, () => g ? tau * randn(rngB) : lapDraw(rngB, b));
      let prev = null, tvSum = 0;
      xs.forEach((x, ix) => {
        const y = evalPoly(bt, x);
        cols[ix].push(y);
        if (prev !== null) tvSum += Math.abs(y - prev);
        prev = y;
      });
      out.tv.push({tv: tvSum, prior});
    }
    cols.forEach((arr, ix) => {
      arr.sort((a, c) => a - c);
      const q = p => arr[Math.min(S - 1, Math.max(0, Math.round(p * (S - 1))))];
      for (const lv of levels) {
        const lo = q(0.5 - lv.p / 2), hi = q(0.5 + lv.p / 2);
        out.yext = Math.max(out.yext, Math.abs(lo), Math.abs(hi));
        out.bands.push({x: xs[ix], lo, hi, lvl: lv.p, prior});
      }
    });
  }
  return out;
}

Code

{
  const clicks = fpDraw;
  const {v, tau, b, deg, xs, levels, bands, yext, lapDraw, evalPoly} = fanData;
  const lam2 = 1 / v;               // Ridge: λ = σ²/τ², σ² = 1
  const lam1 = 2 / b;               // Lasso: λ = 2σ²/b, σ² = 1
  const lines = [];
  for (const prior of ["Gaussian prior (ℓ₂)", "Laplace prior (ℓ₁)"]) {
    const g = prior.startsWith("Gauss");
    const rngL = mulberry32((g ? 51 : 52) + clicks * 7919);
    for (let m = 0; m < 10; m++) {
      const bt = Array.from({length: deg}, () => g ? tau * randn(rngL) : lapDraw(rngL, b));
      for (const x of xs) lines.push({x, y: evalPoly(bt, x), id: prior + "|" + m, prior});
    }
  }
  // paint wide/light first, narrow/dark on top → gradient
  const bandMarks = levels.map(lv =>
    Plot.areaY(bands.filter(d => d.lvl === lv.p),
      {x: "x", y1: "lo", y2: "hi", fx: "prior", fill: lv.c}));
  return html`<div>
    ${themedPlot({height: 320, grid: true,
      fx: {label: null},
      y: {domain: [-1.1 * yext, 1.1 * yext], clamp: true, label: "f(x)"},
      marks: [
        ...bandMarks,
        Plot.ruleY([0], {stroke: "#0f2f5c", strokeDasharray: "3 4", strokeOpacity: 0.6}),
        Plot.line(lines, {x: "x", y: "y", z: "id", fx: "prior",
          stroke: "#c2680c", strokeOpacity: 0.55, strokeWidth: 1})
      ]})}
    <p style="font-variant-numeric: tabular-nums">
      Shared prior variance = ${v.toFixed(3)} &nbsp;→&nbsp; equivalent penalties (σ² = 1):
      &nbsp; <b>Ridge λ = ${lam2.toFixed(2)}</b> &nbsp;|&nbsp; <b>Lasso λ = ${lam1.toFixed(2)}</b>
    </p></div>`;
}

Figure 5: The prior over functions, drawn as a fan: at every x, the nested bands are the symmetric 50/75/90/99% coverage intervals of f(x) under each prior - dark blue core, light blue fringe - with ten sampled functions traced in amber and the two priors matched to the same variance, so whatever differs between the panels is shape of belief, not amount of it. Read the band geometry. The Gaussian fan widens smoothly and proportionally: a hedge in every direction. The Laplace fan, at identical variance, holds a tighter dark core - most of its belief hugs simple, low-magnitude and low-complexity functions - while its pale 99% fringe stretches wider at the edges: rare, sudden complexity. Mostly simple or suddenly wild, little in between. The λ readouts are the penalties these beliefs correspond to (at σ² = 1). You never chose a λ. You chose a worldview.

The plot above asks you to use an eagle eye to see differences. If you can’t see them, let me show you a different way. Compress each draw to a single number - its total variation, \(\int_{-1}^{1} |f'(x)|\,dx\), the total amount of wiggle a function commits over the interval. Simple functions (low magnitude coefficients) have low TV while more complex ones (wigglier with higher magnitude coefficients) have high TV. Look at what each prior believes about that:

Code

{
  const priors = [
    {name: "Gaussian prior (ℓ₂)", color: "#2f64ad"},
    {name: "Laplace prior (ℓ₁)", color: "#c2680c"}
  ];
  const allTv = fanData.tv.map(d => d.tv).sort((a, c) => a - c);
  const xmax = allTv[Math.floor(0.995 * (allTv.length - 1))];
  const grid = linspace(0, xmax, 161);
  const rows = [];
  for (const p of priors) {
    const vals = fanData.tv.filter(d => d.prior === p.name).map(d => d.tv);
    const n = vals.length;
    const mean = vals.reduce((s, x) => s + x, 0) / n;
    const sd = Math.sqrt(vals.reduce((s, x) => s + (x - mean) ** 2, 0) / n);
    const bw = Math.max(1.06 * sd * Math.pow(n, -0.2), xmax / 400);
    for (const x of grid) {
      let dens = 0;
      for (const xv of vals) dens += normPdf(x, xv, bw);
      rows.push({x, y: dens / n, prior: p.name});
    }
  }
  return themedPlot({height: 240, grid: true,
    x: {label: "total variation of f on [−1, 1]", domain: [0, xmax]},
    y: {label: "density"},
    color: {domain: priors.map(p => p.name), range: priors.map(p => p.color), legend: true},
    marks: [
      Plot.areaY(rows, {x: "x", y1: 0, y2: "y", fill: "prior", fillOpacity: 0.18}),
      Plot.line(rows, {x: "x", y: "y", stroke: "prior", strokeWidth: 2})
    ]});
}

Figure 6: The prior’s implied distribution over total variation, one value per sampled function, 8,000 samples per prior at the same matched variance as the fan above. This is the whole disagreement in one picture: the Laplace prior piles its belief hard against low complexity - a tall spike near zero wiggle - then pays for it with a long right tail of occasional wildness. The Gaussian spreads a broad, moderate hump: some wiggle, probably, extreme simplicity and extreme complexity both unlikely. Same variance. Different worldview. Every penalty you have ever tuned was choosing between curves like these.

Sit with the facets for a moment, because this is the reveal underneath the Reveal. The λ you cross-validated was never a number about coefficients. It was a distribution over stories - which shapes of explanation get to be on the table at all. Two penalty families, same mean, same model class, and they believe in visibly different worlds before a single observation lands. Every choice in the modeling pipeline is a statement of this kind, whether its author knows it or not. The rest of this page - and honestly, the rest of the course - is about learning to read those statements.

One more figure, because so far this section has shown only beliefs - and beliefs are supposed to meet data. Here is the full circle: the same energy from the top of the section - likelihood plus regularizer, now over the same degree-9 polynomial family with a Gaussian prior - but this time we watch the posterior over functions as observations accumulate. The truth (amber, dashed) is quadratic; the model doesn’t know that, and it carries ten coefficients of capacity it doesn’t need. Three things to do with your hands. First, push \(N\) up at a sensible prior variance and watch pointwise convergence: the fan tightens onto the true function, everywhere at once. Second, crank the prior variance to its maximum and drop \(N\) below ten - below the coefficient count: with belief effectively switched off and capacity to spare, the posterior mean threads through the data points and the fan detonates between them - interpolation, which you should recognize immediately as mass collapse wearing its function-space costume from Section 3. Third, notice that the variance dial genuinely matters here, unlike in the fan above: it controls how quickly the data can drag the posterior onto the truth - and how badly things go in the meantime when the belief is too weak to carry the load.

Code

viewof pcN = Inputs.range([0.7, 2.7], {value: 1.0, step: 0.05, label: "log₁₀ N (observations)"})
viewof pcTau = Inputs.range([-2, 4], {value: 0, step: 0.1, label: "log₁₀ prior variance τ²"})

Code

{
  const N = Math.round(10 ** pcN);
  const tau2 = 10 ** pcTau;
  // Degree 9 in the raw monomial basis is deliberately chosen for the interpolation
  // demo (10 coefficients → threads ≤10 points when the prior is weak). The monomial
  // Vandermonde is ill-conditioned at d=10; the Cholesky floor should hold, but if
  // the τ²=10⁴ corner renders numeric junk rather than honest detonation, switch the
  // internals to a Legendre basis (narrative unchanged) or cap the τ² slider at 3.
  const sigma = 0.35, deg = 9, d = deg + 1;
  const fTrue = x => 1.5 * x * x - 0.5 * x - 0.2;
  const rng = mulberry32(777);
  const all = Array.from({length: 500}, () => {
    const x = 2 * rng() - 1;
    return {x, y: fTrue(x) + sigma * randn(rng)};
  });
  const data = all.slice(0, N);
  const phi = x => Array.from({length: d}, (_, k) => x ** k);
  // posterior precision: Φ'Φ/σ² + prior⁻¹  (intercept effectively unshrunk: var 1e6)
  const A = Array.from({length: d}, (_, r) =>
    Array.from({length: d}, (_, c) => r === c ? 1 / (r === 0 ? 1e6 : tau2) : 0));
  const bvec = new Array(d).fill(0);
  for (const p of data) {
    const ph = phi(p.x);
    for (let r = 0; r < d; r++) {
      bvec[r] += ph[r] * p.y / sigma ** 2;
      for (let c = 0; c < d; c++) A[r][c] += ph[r] * ph[c] / sigma ** 2;
    }
  }
  const L = chol(A);
  const m = cholSolve(L, bvec);
  const zs = [
    {z: 2.576, c: "#dbeafe"}, {z: 1.645, c: "#a9c9ef"},
    {z: 1.150, c: "#6b9fdc"}, {z: 0.674, c: "#2f64ad"}
  ];
  const xs = linspace(-1, 1, 121);
  const bands = [], meanline = [], truthline = [];
  for (const x of xs) {
    const ph = phi(x);
    const Ainv_p = cholSolve(L, ph);
    const varf = Math.max(ph.reduce((s, v, k) => s + v * Ainv_p[k], 0), 0);
    const mu = ph.reduce((s, v, k) => s + v * m[k], 0);
    const sd = Math.sqrt(varf);
    zs.forEach((zz, iz) => bands.push({x, lo: mu - zz.z * sd, hi: mu + zz.z * sd, iz, c: zz.c}));
    meanline.push({x, y: mu});
    truthline.push({x, y: fTrue(x)});
  }
  const rmse = Math.sqrt(meanline.reduce((s, p, i) => s + (p.y - truthline[i].y) ** 2, 0) / xs.length);
  const bandMarks = zs.map((zz, iz) =>
    Plot.areaY(bands.filter(b => b.iz === iz), {x: "x", y1: "lo", y2: "hi", fill: zz.c}));
  return html`<div>
    ${themedPlot({height: 320, grid: true,
      y: {domain: [-2.5, 3], clamp: true, label: "f(x)"},
      marks: [
        ...bandMarks,
        Plot.line(truthline, {x: "x", y: "y", stroke: "#c2680c", strokeWidth: 2, strokeDasharray: "5 4"}),
        Plot.line(meanline, {x: "x", y: "y", stroke: "#0f2f5c", strokeWidth: 2}),
        Plot.dot(data, {x: "x", y: "y", r: 2.5, fill: "currentColor", fillOpacity: 0.55})
      ]})}
    <p style="font-variant-numeric: tabular-nums">
      N = ${N} &nbsp;|&nbsp; equivalent λ = σ²/τ² = ${(sigma ** 2 / tau2).toExponential(1)}
      &nbsp;|&nbsp; RMSE(posterior mean, truth) = <b>${rmse.toFixed(3)}</b>
    </p></div>`;
}

Figure 7: The posterior over functions for a degree-9 Bayesian polynomial regression: blue bands are the pointwise 50/75/90/99% posterior intervals for f(x) itself (no noise term), navy is the posterior mean, amber-dashed is the quadratic truth, dots are the observations. Sensible prior + growing N: the fan collapses onto the truth pointwise, surplus capacity and all. Prior variance maxed with N below the ten-coefficient count: interpolation - the mean visits every dot and the belief between dots detonates. The convergence is always coming; the prior variance sets the fare, and too weak a prior means paying in wreckage first.

Name what you just did with those two sliders, because the field does it every day under a different vocabulary. Cross-validation is our ritual for exactly this dial: a grid search over λ - over prior variances - hunting for the strength of belief that best balances mass collapse against prior stubbornness. The explanation that is most data-ready. CV finds that balance by loop and holdout, honestly but expensively. Hold the thought, though: Section 8 will show that probability theory can render the same judgment directly, no loop required - and the quantity that does it has been lurking in this section’s corner the whole time.

5. The asymmetry: where priors actually matter

If priors are always present, the next question is when they matter. The answer falls out of arithmetic so simple it’s easy to miss its consequences.

The log-likelihood is a sum over data: \(N\) terms, growing \(O(N)\). The log-prior is one term: \(O(1)\). So as \(N\) grows, the likelihood swamps any fixed prior; the MAP estimate converges to the MLE; your carefully chosen beliefs become a rounding error. Priors are irrelevant, asymptotically - which is to say, priors are irrelevant exactly when you no longer need them.

The formal version is the Bernstein–von Mises theorem (van der Vaart, Asymptotic Statistics, ch. 10). Remember this asymptotic: it is precisely what the “just add data” strategy leans on, a point the book makes much of.

The turn is that modern machine learning does not live in that asymptote. ResNet-50: ~25 million parameters. ImageNet: ~1 million images. We routinely fit systems with far more unknowns than observations - underdetermined by construction, permanently pre-asymptotic. Priors matter precisely where modern ML lives. Feel the tug-of-war yourself:

Figure 8: A single parameter θ, true value 1.2 (dashed). The prior (centered at 0, width τ) fights the likelihood (N observations). Slide N: at small N the posterior hugs the prior; as N grows it gets dragged to the data and the MAP marker converges to the MLE marker. The O(N) sum beats the O(1) term - eventually. Modern ML lives in the ‘not yet.’

6. Real priors: truncate the hypothesis class

Independent Gaussians on weights - the Ridge prior - are the least interesting priors in machine learning. They shrink everything a little and forbid nothing. The priors that do the real work are harsher: hard truncations that assign zero mass to functions that cannot be the answer.

We are very fortunate to live in a world with structure - “peanut butter egg dirt” isn’t a meaningful sentence and a random collection of colored pixels arranged in a grid is not a meaningful image. The world is full of impossibilities and the best models are those that encode that fact.

That is what an architecture is. Convolutional locality says: functions where a pixel’s meaning depends on arbitrary far-away pixels are impossible - deleted, not down-weighted. Weight sharing says: a feature detector that works in the top-left corner works everywhere, and function space shrinks by orders of magnitude at a stroke. Sequential structure, smoothness constraints, bottlenecks - each one is a belief about the world, stated in code instead of in a density, each one deleting a region of hypothesis space outright.

Figure 9: A function space with prior mass as brightness. Flat: mass everywhere, and the data must do all the work. Soft Gaussian: concentrated but nowhere zero - every bad function retains a foothold. Hard truncation: entire regions simply deleted. The volume readout is the point: in high dimensions, deletion is devastatingly efficient - removing half of each of d dimensions removes 1 − 2⁻ᵈ of the volume.

Read the three modes as postures toward belief.

The flat prior is doubt - the refusal to believe anything at all, and Section 3 already showed you where doubt ends: spikes. Data science can only happen here in rare cases.
The soft Gaussian is a hedged bet - probably small, probably smooth, but nothing ever ruled out. Every bad function keeps a foothold it will surrender only under data. But every bad function keeps a foothold it will not surrender until presented with an overwhelming amount of data
Truncation is belief without a hedge: the model may only consider what we hold to be possible, full stop. That sounds like the strictest of the three postures, and it is - but it is also the closest to how the world actually works and humans behave. Nobody believes that their photograph taken with the Wall Street Bull in NYC (you know what photo I’m taking about…) is going to result in a random collection of colors and shapes. Nobody hedges over whether a pixel’s meaning depends on a pixel in a different photograph. You consult your lived experience, declare the combination impossible, and delete it without apology. The question a good architect asks is never “how should I penalize everything a little?” It is: what does my human experience tell me cannot be true here? Truncation isn’t the radical option. It’s the honest transcript of what your brain was already doing.

Take stock of what the last three sections established, because it deserves to be said once, plainly, as the conclusion you just derived rather than a slogan: architectures are priors, losses are priors, training data is a prior. It’s priors all the way down. The only live question - and it is the question of this course - is whose knowledge gets encoded, and with what care.

7. What the point estimate throws away

Everything so far has ended in a point estimate - one \(\hat\theta\), whether MLE or MAP. And a point estimate is a lossy summary of a posterior. What dies in the compression?

Uncertainty: MAP knows the peak’s location and nothing about its width, so a razor-thin posterior and a broad plateau report the same number.
Model comparison: or really model comparison that properly propagates uncertainty under the rules of probability. Concretely: fit two different model classes until each sits at its own peak, both explaining the data about equally well. MAP hands you two numbers and no principled way to rank them - ranking the families requires a quantity that accounts for each class’s entire posterior, not just its best citizen. Hold this thought for a bit.
The one this page has been building toward: the ability to sample: Our stated goal was never just prediction. We want the distribution itself - an object that represents uncertainty, supports probabilistic reasoning, and can be drawn from. A single function is not a generative model. (Second spring: supervised learning can shrug at this and ship the point estimate. Generation cannot - generation is sampling. There is no point-estimate version of “draw me a new dog.”)

Before moving on, see the compression loss with your own eyes. Suppose each density below shows the posterior for the same regression coefficient - say, the effect of a new treatment - arising in three different data situations:

Figure 10: Three posteriors engineered to report the identical MAP - the amber rule at θ̂ = 1 - while being utterly different objects: a razor-thin spike, a broad plateau, and a bimodal density that parks nearly half its mass at a location MAP never mentions. The amber dots along the floor are ten samples from each. A point estimate is a location with no mass, a peak with no shape - and the samples, not the peak, are what generation is made of.

The first panel is the dream: data so informative that the posterior is a spike - given the model, there is next to no uncertainty about the coefficient. The second panel reports the identical MAP and tells an entirely different story: the data are weak, the posterior is broad, and the honest conclusion is “roughly one, but I wouldn’t bet the clinic on it.” MAP cannot tell these two situations apart - and they demand opposite behavior. If θ is a treatment effect, the first posterior licenses decisive action; the second licenses a bigger study. The third panel is the situation that overparameterized models manufacture constantly - multimodality. Same MAP, stranger story: the coefficient is one of two values - strongly positive or moderately negative - and we genuinely do not know which. (Sometimes both modes even yield identical predictions - a lack of uniqueness baked into the estimation problem itself.) Notice that here every single-number summary fails at once. The MAP picks a mode and hides its rival. The posterior mean is worse: it lands between the humps, at a value the model itself considers nearly impossible. No point can carry this story.

And now watch the amber dots as you resample, because this figure is secretly about sampling. Draws from the spike agree with one another - push them through any downstream computation and every simulated future looks the same. Draws from the plateau scatter - propagate them and your predictions inherit the honest wobble. Draws from the bimodal split into two camps - two coherent, incompatible stories about the world, sampled in proportion to how much we believe each. That is what “the posterior is the object” means in practice: anything you do with a fitted model - predict, simulate, decide, generate - happens by drawing from this distribution and pushing the draws through. The samples carry the story. MAP is silent; the dots talk.

So, how do we go from MAP to working with the full posterior? MAP estimates mean that we know (or at least can evaluate the maximum of) the numerator of Bayes’ theorem, \(P(\mathcal D \mid \theta)\, P(\theta)\). Taking the max of this expression tosses information about the posterior shape - we saw this above. What we haven’t dealt with is that MAP throws away information about the only other piece left: the denominator - which has, at last, stepped out of the corner it has been lurking in since Section 4. Ghostface, mask off. We know it by its formal name, the marginal likelihood:

\[ P(\mathcal D) \;=\; \int P(\mathcal D \mid \theta)\, P(\theta)\, d\theta. \]

The marginal likelihood is really the evidence. It doesn’t grade a fitted parameter; it grades a model class: how well does this whole family of hypotheses, weighted by its prior, explain the data? How much posterior mass exists at all? How much fine-tuning did the fit require? MAP proportionality let us cancel it, so we did, for three sections. That free ride is over.

Note how this concept links back to the idea of “priors on functions” - our prior establishes the set/continuum of functions that are possible and assigns a “before conditioning on data” probability to each one based on our understanding of the world. We then condition our prior beliefs on the data that we observed - algorithmically, via Bayes’ theorem, which is always true. But our prior specification put some mass on each possibility meaning that the data we observed could have plausibly been generated by any of those functions. The marginal likelihood, then, is the principled way to turn this distribution over functions into a single number that grades the model class as a whole - taking our prior beliefs into account, what is the prior weighted likelihood (read as expected value) that we would observe the data we actually saw? Written as the expectation it secretly is:

\[ P(\mathcal D) \;=\; \mathbb E_{\theta \sim P(\theta)}\big[\, P(\mathcal D \mid \theta) \,\big]. \]

In plain English - the marginal likelihood quantifies how surprised or unsurprised we would be to see the data we did given our prior beliefs (both encoded via the explicit prior term and encoded through the choice of architecture).¹

8. The probability pivot

Seems great, right? This quantity is actually more familiar than you probably realize! If we alter the framework around the marginal likelihood just a little to shift the uncertainty from the function itself (“What is the probability that any one of these functions generated the data?”) to uncertainty about what the training data would look like if I drew a sample of size \(N\) a second time (“What is the probability that this function would generate a slightly different training set?”), we get our old friend cross-validation. Let’s establish this near equivalence more formally.

Some notation first. Let \(\mathbf y\) be the outcome. In the supervised case, this can be numbers, class labels, vectors, matrices, or whatever else your model predicts. In the unsupervised case, \(\mathbf y\) is just a collection of data points. Let \(\mathbf x\) be a corresponding set of features (a.k.a. things to be conditioned on and assumed known). Let \(\boldsymbol \theta\) be a collection of model parameters that are going to be learned from the data - think of these as regression coefficients, model weights, etc. that are learned via an optimization procedure. Finally, let \(\boldsymbol \lambda\) be a collection of hyperparameters that must be specified a priori in order for the estimation/optimization procedure to learn values of \(\boldsymbol \theta\). These are parameters that are unknown but must be specified in order for the optimization procedure to be tractable - think the regularization parameter in Ridge regression and LASSO. But, this set of hyperparameters also includes the choice of architecture itself - the type of functional form to be learned (linear models, step functions, splines, piecewise linear compositions, etc.) as well as any parameters required to encode the prior (what shape does it take?, what values are impossible?, etc.).

When we see our training data, we see all \(N\) points at once. But, let’s pretend it didn’t happen this way. Instead, imagine that your data arrived as a stream - we saw \(y_1\) before \(y_2\) before \(y_3\) and so on. When each new data point arrives, we use any previously seen data to try to predict the new outcome (using only the prior to predict \(y_1\) since we didn’t see anything before that one). Predict \(y_1\) before seeing anything, using the prior alone. Given our prior beliefs on \(\boldsymbol \theta\), we can define a full predictive distribution for a new outcome conditioning on everything seen before as:²

\[ P(y^*_i \mid y_{<i}, \mathbf x_{\le i}, \boldsymbol \lambda) \;=\; \int P(y^*_i \mid \mathbf x^*_i, \boldsymbol \theta, \boldsymbol \lambda)\, P(\boldsymbol \theta \mid y_{<i}, \mathbf x_{<i}, \boldsymbol \lambda)\, d\boldsymbol \theta . \]

Then, plugging in the actual \(y_i\) and evaluating its likelihood from this distribution gives us a predictive likelihood. Note that once we plug in the observed \(y_i\), this value is a constant - a number, not a distribution. The predictive distribution was the machine; the predictive likelihood is its output.

In words, the predictive likelihood for \(y_i\) evaluated in the stream is the likelihood we would observe \(y_i\) given some values on the model parameters \(\boldsymbol \theta\) (the only uncertain part under this construction), the model hyperparameters \(\boldsymbol \lambda\), and the values of any conditioners \(\mathbf x_i\) weighted by the posterior distribution on \(\boldsymbol \theta\) given all previously observed data points and conditioners. We’re making sequential predictions (nothing new) but we also account for the uncertainty in our model parameters by marginalizing over the posterior distribution of \(\boldsymbol \theta\). Marginalizing over beliefs is the key tenet of Bayesian inference that allows for proper handling of uncertainty!

Note that this is not all that different from what we do when we use the \(K\)-fold cross validation estimator to estimate out-of-sample predictive performance - our data arrives in two chunks, the training set and the validation set, and we try to predict the validation set given our whole set of features and \(\mathbf y\) from the training chunk.³ This process just takes that idea to the extreme!

Using this recipe, we can define a score for the sequential predictive likelihood of the data stream as the sum of log predictive likelihoods (like LOOCV):

\[ \sum_{i=1}^N \log P(y_i \mid y_{<i}, \mathbf x_{\le i}, \boldsymbol \lambda) \]

Code

{
  const i = pvStep;
  const seen = pivotData.slice(0, i), next = pivotData[i];
  // running ledgers
  let ledgerA = 0, ledgerB = 0;
  for (let k = 0; k < i; k++) {
    const pa = predictive(pivotData, k, 1, 2.0, pivotData[k].x);
    const pb = predictive(pivotData, k, 6, 8.0, pivotData[k].x);
    ledgerA += logNormPdf(pivotData[k].y, pa.mean, pa.sd);
    ledgerB += logNormPdf(pivotData[k].y, pb.mean, pb.sd);
  }
  const grid = linspace(-2.2, 2.2, 88);
  const bandA = grid.map(x => { const p = predictive(pivotData, i, 1, 2.0, x); return {x, lo: p.mean - 2 * p.sd, hi: p.mean + 2 * p.sd, mid: p.mean, which: "Line"}; });
  const bandB = grid.map(x => { const p = predictive(pivotData, i, 6, 8.0, x); return {x, lo: p.mean - 2 * p.sd, hi: p.mean + 2 * p.sd, mid: p.mean, which: "Flexo"}; });
  const marks = [
    Plot.areaY(bandA, {x: "x", y1: "lo", y2: "hi", fill: theme.blue, fillOpacity: 0.12}),
    Plot.line(bandA, {x: "x", y: "mid", strokeWidth: 2, stroke: theme.blue}),
    Plot.areaY(bandB, {x: "x", y1: "lo", y2: "hi", fill: theme.amber, fillOpacity: 0.10}),
    Plot.line(bandB, {x: "x", y: "mid", strokeWidth: 2, stroke: theme.amber, strokeDasharray: "5 3"}),
    Plot.dot(seen, {x: "x", y: "y", r: 3.5, fill: theme.ink, fillOpacity: 0.8})
  ];
  if (next) marks.push(Plot.dot([next], {x: "x", y: "y", r: 7, symbol: "star", fill: theme.amber, stroke: "white", strokeWidth: 1.5}));
  return html`<div>
    ${themedPlot({height: 300, grid: true, y: {domain: [-5.5, 5.5]}, marks})}
    <p style="font-variant-numeric: tabular-nums">
      Sequential score after ${i} points - <b>Line: ${ledgerA.toFixed(2)}</b> &nbsp; <b>Flexo (dashed): ${ledgerB.toFixed(2)}</b>
      ${next ? html`<br><small>★ = the next point, about to be predicted by both.</small>` : html`<br><small>Stream complete.</small>`}
    </p></div>`;
}

Figure 11: Data arrive one at a time. Before each arrival, two models state a full predictive distribution: honest Line (linear, sensible prior) and Flexo (degree-6 polynomial, diffuse prior - the kind of model that aces training sets). The ledger is the running sum of log predictive scores. Step to the end: Flexo fits every point it has seen beautifully - and loses, because prediction is scored before seeing, and confident wrongness compounds. Occam’s razor, emerging from bookkeeping, with no penalty term anywhere in sight.

This sequential predictive score works much like a regularization penalty: the training data may fit the model class well, but each new data point that comes in needs to also be explained by the same model class!

That still looks like cross-validation with extra steps, right? And you’d be right - until probability theory collects its rent. Write the score down and stare at it:

\[ \sum_{i=1}^{N} \log P(y_i \mid y_{<i}, \mathbf x_{\le i} , \boldsymbol \lambda) \;=\; \log \prod_{i=1}^{N} P(y_i \mid y_{<i}, \mathbf x_{\le i}, \boldsymbol \lambda) \;=\; \log P( y_1 , y_2, ... , y_N \mid \mathbf x, \boldsymbol \lambda) \]

The chain rule collapses the whole sequential game into a single joint probability - the probability that we would see the data we did given conditioners and a fixed model class \(\boldsymbol \lambda\) with the unknown parameters marginalized out. This is a quantity that is essentially cross validation. A quantity that has an explicit analytical formula that does not require \(K\) different training runs. A quantity that inherently penalizes against models that overfit or underfit the data.

That joint probability, with the parameters integrated out at every step, is the marginal likelihood. The denominator of the Bayes’ Theorem formula and the honest one-step-ahead scorekeeping you just watched are the same number. The evidence is leave-future-out cross-validation, computed by probability theory instead of by loop. That is what it’s worth: the marginal likelihood is not a normalizing nuisance. It is the compass used throughout machine learning!

And look back at the widget’s ledger. The running number you watched accumulate for Line and Flexo was the log marginal likelihood of each model class the entire time - you were not watching a demonstration of scorekeeping, you were watching two evidences being computed. Which also answers the thought Section 7 asked you to hold: that is how you rank model families - not by comparing peaks, but by comparing how well each class, weighted over its whole posterior, predicted the stream.

The identity is in the literature: Fong & Holmes, “On the marginal likelihood and cross-validation” (Biometrika 2020). For evidence-as-Occam: MacKay ITILA ch. 28; Kass & Raftery, “Bayes Factors” (JASA 1995).

9. Generative ML is Bayes with the training wheels off

The marginal likelihood is the golden goose for modern ML - a single equation that captures the ability of your model to generalize.⁴ So this final act runs on three questions. What, in full generality, is the optimization that machine learning is actually trying to do? Why does the field mostly refuse to compute the very quantity it agrees on? And why does generative modeling - the subject this whole page has been winding toward - get no such refusal?

One more thing before the bad news. Look one last time at the shape of the machine Section 8 built: predict the next thing, given everything seen so far, and score the whole run by summed log predictive likelihood. Remember that shape. It is about to become the most valuable object in machine learning.

9.1 The inner and outer problems

Start from the practitioner’s honest situation. You have a million values your hyperparameters could take and a million architectures you could try - and that is before anyone mentions the prior. The working paradigm of the field, empirical risk minimization, handles half of this cleanly: given a setting of \(\boldsymbol \lambda\), the parameters are solved as an argmin,

\[ \hat{\boldsymbol \theta}(\boldsymbol \lambda) \;=\; \arg\min_{\boldsymbol \theta} \; \Big[ \text{loss}(\boldsymbol \theta ; \mathcal D) + \text{regularizer}(\boldsymbol \theta ; \boldsymbol \lambda) \Big]. \]

Call this the inner problem - it is conditioned on \(\boldsymbol \lambda\) from the start, and it is the part we are good at. But the conditioning is exactly the issue: someone still has to set those other values. How do we choose between architectures? How do we set the tuning parameters? How do we pick the prior’s shape?

The principled answer assembles itself from everything this page has built: take the inner problem you conditioned on \(\boldsymbol \lambda\) and wrap it in the marginal likelihood’s integral so the parameters are handled by marginalization rather than by faith in a single argmin, and maximize the result with respect to \(\boldsymbol \lambda\):

\[ \hat{\boldsymbol \lambda} \;=\; \arg\max_{\boldsymbol \lambda} \; \log \int \exp\Big\{ -\Big[ \text{loss}(\boldsymbol \theta ; \mathcal D) + \text{regularizer}(\boldsymbol \theta ; \boldsymbol \lambda) \Big] \Big\}\, d\boldsymbol \theta. \]

The integrand is not decoration. By Section 4’s energy identity, \(\exp\{-[\text{loss} + \text{regularizer}]\} \propto P(\mathcal D \mid \boldsymbol \theta, \boldsymbol \lambda)\, P(\boldsymbol \theta \mid \boldsymbol \lambda)\) - so the thing being integrated is exactly likelihood times prior, and the outer objective is exactly the log marginal likelihood. The inner problem survives inside it as the integrand’s peak: once \(\hat{\boldsymbol \lambda}\) is chosen, \(\hat{\boldsymbol \theta}(\hat{\boldsymbol \lambda})\) is recovered by the familiar argmin. An inner problem for parameters, an outer problem for everything else, with the evidence as the bridge between them.

This two-level formulation is, more or less, machine learning as a whole. And it is remarkably uncontroversial - frequentist or Bayesian, theorist or engineer, nearly everybody agrees that this is the problem. So why doesn’t everybody simply compute the outer objective and go home?

9.2 The integral is hard - and where cross-validation goes to die

Believe me when I say this because I’ve spent a long time working on this:

This. Integral. Is. Hard. To. Compute.

The evidence \(\int P(\mathcal D\mid\boldsymbol\theta,\boldsymbol\lambda)\,P(\boldsymbol\theta\mid\boldsymbol\lambda)\,d\boldsymbol\theta\) has closed forms in a museum’s worth of conjugate special cases and essentially nowhere else. In high dimensions, quadrature dies, Monte Carlo estimates arrive with variances you could park a truck in, and the integrand is a needle of posterior mass threading a haystack of prior volume. The most important quantity in the course is, in general, unobtainable.

So what do we do? For supervised methods, we can go back to cross-validation. This is the entire reason \(K\)-fold still rules supervised practice - not superiority, availability. The loop is possible when the integral is not: refitting a model \(K\) times asks for nothing but compute and patience, and Section 8 showed the two quantities are chasing the same target anyway. A reasonable trade, in the regime where it is a trade at all.

But there are two scenarios where the loop stops making sense - and they happen to be the two scenarios that define modern ML.

Supervised learning at scale. The marginal likelihood would be genuinely useful for the largest models we build - and it is precisely there that its stand-in becomes absurd. Training runs for frontier systems - the trillion-parameter class, Claude Fable 5 and its peers - are measured in months of wall-clock time across tens of thousands of accelerators, with price tags that get reported with the word “billion” nearby. Evaluating a validation fold is cheap. Now propose, with a straight face, retraining the model from scratch \(K-1\) more times to rotate the folds. The loop’s only virtue was availability, and at scale, availability is gone too.
No clean answer to “was this generation successful?” Say we trained a generator and drew a sample. How do we judge whether that sample succeeded? How do we even split data into folds - and what would we do with the validation fold once we had it? There is no clear right and wrong for a generated dog - at least nothing beyond human eye tests and other ad hoc rituals - so what, exactly, is being measured when we “cross-validate” a generator? The loop does not merely get expensive here. Its question stops parsing.

So we need other solutions - good approximations to the marginal likelihood, which is exactly what the escapes of Section 9.4 are. But for generative methods there is a better reason still, and it is the reason this page exists: the marginal likelihood is the quantity that exactly answers the question we are asking when we build a generative model. Not a proxy for it. The thing itself. We will make that claim precise in a moment.

And notice, before we do, that the three springs this page wound have just released. Collapse (Section 3): with no prediction task, nothing papers over a model that put its mass in the wrong place - samples come from wherever the mass is. Uncertainty (Section 7): there is no point-estimate version of a new draw - generation is sampling. And the evidence: supervised learning could cancel the denominator and loop instead; generative modeling just watched both escape routes close. Generative ML is the fundamental problem of Section 2 with no place left to hide.

9.3 What a generative model formally is

By now you have all the parts, so the formal statement will read less like a definition and more like a confession of what you already believe.

We want a random draw of \(\mathbf x\) from \(f(\mathbf x)\) - the big ol’ distribution in the sky that generates the thing that we see. Maybe it’s pictures of dogs. Maybe it’s videos of farting wombats. Maybe it’s a story about a gruff southern father finally telling his son that he loves him. Different things, but all assumed to be generated from some unknown joint distribution over pixels/frames/words.

Three facts about \(f\):

We do not know it (too complicated).
It is impossible to write down by hand (even more complicated still).
It is a proper probability density: nonnegative, integrating to one, its mass distributed over a space too large to think about directly. We can assume this to be true because we are the creators.

That is the entire problem statement of modern generative AI.

Make it concrete twice.

Images: \(f\) is a density over pixel space - a 256×256 RGB image is one draw from a distribution on \(\mathbb R^{196{,}608}\). Try sampling that space with no structural knowledge:

Code

{
  npClicks; // reactive dependency
  const w = 160, h = 160;
  const canvas = DOM.canvas(w, h);
  const ctx = canvas.getContext("2d");
  const img = ctx.createImageData(w, h);
  for (let k = 0; k < img.data.length; k += 4) {
    img.data[k] = 256 * Math.random(); img.data[k+1] = 256 * Math.random();
    img.data[k+2] = 256 * Math.random(); img.data[k+3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  return html`<div style="display:flex; gap:12px; align-items:center;">
    ${canvas}
    <div style="max-width:340px;font-size:13px;line-height:1.4;opacity:.75">
      Every photograph you have ever cared about is a draw from <em>f</em> - and not one of them looks like the panel on the left.  That gap is the whole problem.
    </div></div>`;
}

Figure 12: Uniform draws from pixel space. Press it as many times as you like: static, every time. Nearly all of pixel space is noise; every image anyone has ever cared about lives on a vanishing wisp of structure threading through it. f puts essentially all its mass on the wisp - and we cannot write down where the wisp is. (When we build discrete image generators, a button just like this one will make the same point in token space.)

Text: \(f\) is a distribution over word sequences - every sentence ever written, one draw each, weighted by how the language actually goes:

Code

{
  const steps = [
    {prefix: "The …", bars: [["dog", .18], ["sun", .09], ["problem", .07], ["quick", .04], ["…thousands more", .62]]},
    {prefix: "The dog …", bars: [["chased", .22], ["barked", .18], ["slept", .12], ["ran", .11], ["…", .37]]},
    {prefix: "The dog chased the …", bars: [["squirrel", .34], ["ball", .22], ["mailman", .13], ["cat", .12], ["the", .002], ["…", .188]]},
    {prefix: "The dog chased the squirrel …", bars: [["up", .3], ["across", .2], ["into", .17], ["yesterday", .05], ["…", .28]]}
  ];
  const s = steps[nwStep];
  return html`<div>
    <p style="font-size:1.05em"><b>${s.prefix}</b></p>
    ${themedPlot({height: 190, marginLeft: 90, x: {label: "P(next word | so far)"},
      marks: [Plot.barX(s.bars.map(([w, p]) => ({w, p})), {y: "w", x: "p", sort: {y: "-x"}, fill: theme.blue, rx: 2})]})}
  </div>`;
}

Figure 13: A distribution over sentences, expressed one conditional at a time. (Probabilities hand-authored for illustration - no model behind this, which is rather the point: you can feel the shape of f without being able to write it down. Note the form: next word, given the words so far. You have seen that form before, one section ago.)

What we are actually trying to learn is the underlying joint distribution of the features themselves - \(f\) as an object. What we get are samples from it. Potentially billions or trillions of samples - but samples all the same. And Section 3 told you, at slider-length, exactly what happens if we fit a distribution to samples without imposing prior constraints: mass collapse. Delta spikes on the training set, at production scale. Even with millions upon millions of spikes, our goal is still to fill in the area between those spikes - the empirical distribution alone can only ever re-emit what it has already seen. Billions of observed dogs, and zero new ones. Nothing generative about that.

But, Kevin, what if I have millions of examples of dog images? Won’t the empirical distribution converge to the true distribution? Yes, it will - eventually, and “eventually” is doing back-breaking work in that sentence. The number of examples needed to fill the space of all possible images of dogs is astronomical. The number of possible images is \(256^{256 \times 256 \times 3}\), and the number of images of dogs is a vanishingly small fraction of that. In fact, Stone’s minimax theory (Stone 1982) tells us that the sample size needed to estimate a density to fixed accuracy grows exponentially in the dimension of the space - the curse of dimensionality, stated as a theorem. Even a billion dog images samples an infinitesimal fraction of pixel space.

So we place priors - and in the generative case, the priors finally make vivid sense. To specify a distribution over images of dogs, we have to hand the model enough of the world to say: the only distributions on the table are ones where there is an object with eyes, fur, and a snout - two ears, left and right, above the eyes. And we need to enforce some sense of smoothness: the numerical representation of an image of a chihuahua needs to be close to other representations of chihuahuas. Decline to impose that structure and there are exactly two roads. Either the model must learn those nuances from scratch - which demands enormous flexibility and massive sample sizes, and is, in one sentence, the attention-based approach - or the model floats hopelessly through ultra-high-dimensional pixel space until it settles on some saddle point that emits random images with no respect for the laws of the universe. You pressed the button that samples from that second road a few paragraphs ago.

More formally, the object that grades a trained generator on a single observation is

\[ f(\mathbf x \mid \boldsymbol \lambda) \;=\; \int f(\mathbf x \mid \boldsymbol \theta, \boldsymbol \lambda)\, f(\boldsymbol \theta \mid \boldsymbol \lambda)\, d\boldsymbol \theta \]

the probability that we would observe this sample image of a dog, given our model class, accounting for uncertainty in the learned parameters. That accounting is not a nicety. It is what makes the model bleed each training image into its neighborhood of pixel space, inducing the locality a generator needs - if a picture of a shih tzu lives here, then a little to the left should also be a picture of a shih tzu.

Read the equation in plain English: given the prior knowledge of how the world works that we baked into the model, and what we learned from the data, how likely is it that our distribution generated the things we observed? You know this quantity. It is the marginal likelihood - and it answers the generative question exactly. No stand-in, no side-step, no “something close.” It is not a proxy for the thing we are trying to measure. It is the thing we are trying to measure.

The maximization problem, stated formally. Find the model class \(\boldsymbol \lambda\) - the architecture, the prior, the hyperparameters, exactly the object Section 8’s score grades - that maximizes the marginal likelihood of the observed data:

\[ \boldsymbol \lambda^\star \;=\; \arg\max_{\boldsymbol \lambda} \; \log \int P(\mathcal D \mid \boldsymbol \theta, \boldsymbol \lambda)\, P(\boldsymbol \theta \mid \boldsymbol \lambda)\, d\boldsymbol \theta, \]

with \(\boldsymbol \theta\) standing for everything unknown and integrated over (parameters, latent variables - the structure varies, the shape of the problem doesn’t). Outside conjugate toys, this has no closed form. We cannot evaluate the objective we most want to climb, let alone differentiate our way up it.

The modern field’s answer is a set of computational proxies for exactly this maximization - the ELBO, diffusion, autoregression - and they are the next section. But before any machinery, this quantity is what tells us how successful our generative model is - using the trained model, how likely are the data points that we observed under this model (accounting for our prior beliefs)? Ergo, our goal is to come up with a model that maximizes the marginal likelihood given the training data. Modern generative AI that produces incrementally more humanlike text - the new model increased the marginal likelihood. An image generator that produces higher fidelity images of dogs wearing watermelon hats - the new model increased the marginal likelihood. If the marginal likelihood is hard to compute, then how do we know how to develop models that increase it? The rules of probability hand us two honest levers that give us a route to increase the marginal likelihood associated with a particular model:

Lever 1: increase \(N\). With a sufficiently flexible approximator, more data almost surely increases the likelihood component of the evidence. “Sufficiently flexible” sounds like a demanding hypothesis; it is table stakes. Modern ML runs entirely on universal approximators - \(k\)-nearest neighbors, trees and their ensembles, neural networks - families rich enough to represent essentially anything given enough data. The field cleared this bar decades ago. What increasing \(N\) costs is a different question.

Lever 2: decrease wasted prior volume. Every unit of prior mass parked where the answer isn’t is likelihood the evidence integral can never collect. Sharpen the prior toward the truth - better yet, delete regions the truth cannot occupy - and the same data buys a bigger number.

Don’t take my word for it. Let’s actually build a generative model - the smallest one that still deserves the name - and watch both levers work on it.

The setup: Hidden in a blank unit square is a secret animal - one of five (cat, dog, horse, frog, monkey), possibly mirrored, chosen at random and always contained in the top right corner of the grid (what exactly that means is up to you). The true distribution is uniform over that animal’s silhouette; a “generated animal” is the shape a few hundred points make when drawn from a density, which is not a costume at all: a point cloud genuinely is a sample from a distribution, and a distribution is the thing we are trying to learn. The model chops the square into a \(32 \times 32\) grid and calls a point’s cell a draw from a multinomial - the histogram version of density estimation, no smoothness hiding anywhere - with a Dirichlet prior over the cell probabilities to ensure that the marginal likelihood is defined and computable.

Three types of priors:

The flat prior: concentration spread evenly, the animal could be anywhere.
A soft prior you shape - a locked Gaussian bump holds the center (“it’s centered” is not optional), and you may add up to five more at a width you control; every cell keeps positive mass, so this prior can never be falsified, only out-predicted.
A hard mask, carved in coarse \(4 \times 4\) blocks: it starts as the upper-right \(3 \times 3\) of those blocks, and you add or delete blocks according to your understanding of the phrase “top right” - every deleted block is a declaration of impossible.

All three become non-symmetric Dirichlets whose per-cell concentration tracks the belief you drew, all three are scored on the same observed points, and all three evidences are exact and closed-form. Play with this widget a bit and think about what you see in terms of the marginal likelihood and its relationship with sample size and generation quality.

Code

// "Draw me an animal" - three priors (flat / soft MVN mixture / hard mask), exact
// Dirichlet-multinomial evidence, shared observed points. Fully client-side, zero assets.
daShapes = {
  const inCircle = (x, y, cx, cy, r) => (x - cx) ** 2 + (y - cy) ** 2 <= r * r;
  const inEll = (x, y, cx, cy, rx, ry) => ((x - cx) / rx) ** 2 + ((y - cy) / ry) ** 2 <= 1;
  const inRect = (x, y, x0, x1, y0, y1) => x >= x0 && x <= x1 && y >= y0 && y <= y1;
  const inTri = (px, py, ax, ay, bx, by, cx, cy) => {
    const s1 = (bx - ax) * (py - ay) - (by - ay) * (px - ax);
    const s2 = (cx - bx) * (py - by) - (cy - by) * (px - bx);
    const s3 = (ax - cx) * (py - cy) - (ay - cy) * (px - cx);
    return (s1 >= 0 && s2 >= 0 && s3 >= 0) || (s1 <= 0 && s2 <= 0 && s3 <= 0);
  };
  // all silhouettes live inside [0.17, 0.83]^2 (the middle 16x16 of a 24-grid)
  const dog = (x, y) =>
    inEll(x, y, .45, .45, .20, .12) || inCircle(x, y, .70, .60, .085) ||
    inEll(x, y, .80, .58, .055, .04) || inTri(x, y, .62, .65, .65, .77, .69, .65) ||
    inRect(x, y, .31, .35, .22, .38) || inRect(x, y, .41, .45, .22, .38) ||
    inRect(x, y, .51, .55, .22, .38) || inRect(x, y, .59, .63, .22, .38) ||
    inCircle(x, y, .255, .52, .038);
  const cat = (x, y) =>
    inEll(x, y, .47, .40, .14, .155) || inCircle(x, y, .585, .62, .092) ||
    inTri(x, y, .52, .68, .55, .81, .585, .69) || inTri(x, y, .60, .69, .635, .81, .665, .68) ||
    inCircle(x, y, .335, .45, .033) || inCircle(x, y, .31, .53, .03) || inCircle(x, y, .30, .61, .028);
  const horse = (x, y) =>
    inEll(x, y, .47, .52, .19, .10) || inEll(x, y, .67, .63, .065, .092) ||
    inEll(x, y, .72, .72, .07, .045) || inTri(x, y, .69, .76, .71, .82, .73, .76) ||
    inRect(x, y, .32, .35, .26, .44) || inRect(x, y, .41, .44, .26, .44) ||
    inRect(x, y, .51, .54, .26, .44) || inRect(x, y, .59, .62, .26, .44) ||
    inEll(x, y, .275, .55, .032, .08);
  const frog = (x, y) =>
    inEll(x, y, .50, .42, .175, .12) || inEll(x, y, .50, .54, .12, .075) ||
    inCircle(x, y, .43, .61, .042) || inCircle(x, y, .57, .61, .042) ||
    inEll(x, y, .33, .37, .07, .092) || inEll(x, y, .67, .37, .07, .092) ||
    inRect(x, y, .27, .40, .28, .32) || inRect(x, y, .60, .73, .28, .32);
  const monkey = (x, y) =>
    inEll(x, y, .49, .44, .12, .15) || inCircle(x, y, .49, .64, .088) ||
    inCircle(x, y, .42, .70, .032) || inCircle(x, y, .56, .70, .032) ||
    inRect(x, y, .365, .40, .28, .44) || inRect(x, y, .58, .615, .28, .44) ||
    inCircle(x, y, .63, .42, .032) || inCircle(x, y, .66, .49, .03) || inCircle(x, y, .67, .57, .028);
  return {list: [["dog", dog], ["cat", cat], ["horse", horse], ["frog", frog], ["monkey", monkey]]};
}

Code

// the secret: species + mirror + a big reservoir of truth points (reseeded by the button)
daSecret = {
  daNewBtn;
  const G = daG, Nres = 100000;
  const pick = mulberry32(9000 + (daNewBtn || 0) * 77 + 13);
  const idx = Math.floor(pick() * daShapes.list.length);
  const mirror = pick() < 0.5;
  const [name, raw] = daShapes.list[idx];
  // shapes are authored in ~[0.17,0.83]^2; place them in the top-right quadrant [0.5,1]x[0.5,1].
  // map quadrant point -> native space: u = (x-0.5)/0.5, then into [0.17,0.83]
  const toNative = (x, y) => ({ u: 0.17 + ((x - 0.5) / 0.5) * 0.66, v: 0.17 + ((y - 0.5) / 0.5) * 0.66 });
  const shape = (x, y) => {
    if (x < 0.5 || y < 0.5) return false;
    const {u, v} = toNative(x, y);
    return mirror ? raw(1 - u + 0.34, v) : raw(u, v);   // mirror about the quadrant's own center
  };
  const rng = mulberry32(1789 + idx * 131 + (mirror ? 7 : 0));
  const pts = [];
  while (pts.length < Nres) { const x = rng(), y = rng(); if (shape(x, y)) pts.push({x, y}); }
  const cellOf = (x, y) => Math.min(G - 1, Math.floor(y * G)) * G + Math.min(G - 1, Math.floor(x * G));
  const ptCell = pts.map(p => cellOf(p.x, p.y));
  const ghost = [];
  for (let r = 0; r < 48; r++) for (let c = 0; c < 48; c++) {
    const x = (c + .5) / 48, y = (r + .5) / 48;
    if (shape(x, y)) ghost.push({x1: c / 48, x2: (c + 1) / 48, y1: r / 48, y2: (r + 1) / 48});
  }
  return {G, name, mirror, pts, ptCell, ghost};
}

Code

// shared cell-grid click editor factory. mode "modes": store click centers (cap 6).
// mode "mask": toggle cells, start ALL ON.
daEditor = function(kind, caption) {
  const G = daG, size = 220, cs = size / G;
  const B = 4, SG = G / B;   // 8x8 super-grid; each super-cell is a 4x4 block
  const setBlock = (a, sc, sr, v) => { for (let r = sr * B; r < sr * B + B; r++) for (let c = sc * B; c < sc * B + B; c++) a[r * G + c] = v; };
  const state = kind === "mask" ? (() => { const a = new Uint8Array(G * G); for (let sr = SG - 3; sr < SG; sr++) for (let sc = SG - 3; sc < SG; sc++) setBlock(a, sc, sr, 1); return a; })() : [{c: G / 2, r: G / 2, locked: true}];
  const canvas = html`<canvas width=${size} height=${size} style="border:1px solid currentColor;border-radius:4px;cursor:crosshair;touch-action:none">`;
  const ctx = canvas.getContext("2d");
  const root = html`<div><div style="font-size:.78em;font-variant-caps:all-small-caps;letter-spacing:.04em;opacity:.65;margin:0 0 3px 4px">${caption}</div></div>`;
  const rc = i => ({c: i % G, r: Math.floor(i / G)});
  const draw = () => {
    ctx.clearRect(0, 0, size, size);
    // top-right 16x16 guide (where animals live on the 32-grid)
    ctx.fillStyle = "rgba(128,128,128,0.06)";
    ctx.fillRect(16 * cs, 0, 16 * cs, 16 * cs);
    if (kind === "mask") {
      ctx.fillStyle = "rgba(194,104,12,0.26)";
      for (let sr = 0; sr < SG; sr++) for (let sc = 0; sc < SG; sc++) if (state[(sr * B) * G + sc * B]) ctx.fillRect(sc * B * cs, (G - (sr + 1) * B) * cs, B * cs, B * cs);
      ctx.strokeStyle = "rgba(128,128,128,0.4)";
      for (let k = 0; k <= SG; k++) { ctx.beginPath(); ctx.moveTo(k * B * cs, 0); ctx.lineTo(k * B * cs, size); ctx.stroke(); ctx.beginPath(); ctx.moveTo(0, k * B * cs); ctx.lineTo(size, k * B * cs); ctx.stroke(); }
    } else {
      for (const m of state) {
        const px = (m.c + .5) * cs, py = (G - 1 - m.r + .5) * cs;
        if (m.locked) { ctx.strokeStyle = "rgba(47,100,173,0.95)"; ctx.lineWidth = 2; ctx.beginPath(); ctx.arc(px, py, cs * .42, 0, 7); ctx.stroke(); }
        else { ctx.fillStyle = "rgba(47,100,173,0.9)"; ctx.beginPath(); ctx.arc(px, py, cs * .42, 0, 7); ctx.fill(); }
      }
    }
    ctx.strokeStyle = "rgba(128,128,128,0.22)";
    for (let k = 0; k <= G; k++) { ctx.beginPath(); ctx.moveTo(k * cs, 0); ctx.lineTo(k * cs, size); ctx.stroke(); ctx.beginPath(); ctx.moveTo(0, k * cs); ctx.lineTo(size, k * cs); ctx.stroke(); }
  };
  const emit = () => { root.value = kind === "mask" ? Array.from(state, v => !!v) : state.map(m => ({c: m.c, r: m.r})); root.dispatchEvent(new Event("input", {bubbles: true})); };
  const at = ev => { const b = canvas.getBoundingClientRect(); const c = Math.max(0, Math.min(G - 1, Math.floor((ev.clientX - b.left) / cs))); const rr = Math.max(0, Math.min(G - 1, Math.floor((ev.clientY - b.top) / cs))); return {c, r: G - 1 - rr, i: (G - 1 - rr) * G + c}; };
  let painting = false, pv = 0;
  canvas.addEventListener("pointerdown", ev => {
    const {c, r, i} = at(ev);
    if (kind === "mask") { painting = true; canvas.setPointerCapture(ev.pointerId); const sc = Math.floor(c / B), sr = Math.floor(r / B); pv = state[(sr * B) * G + sc * B] ? 0 : 1; setBlock(state, sc, sr, pv); draw(); }
    else {
      const hit = state.findIndex(m => m.c === c && m.r === r);
      if (hit >= 0) { if (!state[hit].locked) state.splice(hit, 1); }
      else if (state.length < 6) state.push({c, r});
      draw(); emit();
    }
  });
  canvas.addEventListener("pointermove", ev => { if (!painting) return; const {c, r} = at(ev); const sc = Math.floor(c / B), sr = Math.floor(r / B); if (state[(sr * B) * G + sc * B] !== pv) { setBlock(state, sc, sr, pv); draw(); } });
  canvas.addEventListener("pointerup", () => painting = false);
  root.append(canvas);
  if (kind === "mask") {
    const kept = html`<span style="font-size:.78em;opacity:.7"></span>`;
    const saveBtn = html`<button style="margin-top:4px">Save mask &amp; recompute</button>`;
    const refreshKept = () => kept.textContent = ` ${state.reduce((s, v) => s + v, 0)} cells kept`;
    canvas.addEventListener("pointerup", refreshKept);
    canvas.addEventListener("pointerdown", refreshKept);
    saveBtn.onclick = () => { emit(); refreshKept(); };
    root.append(html`<div style="display:flex;align-items:center;gap:8px">${saveBtn}${kept}</div>`);
    refreshKept();
  }
  draw(); emit();
  return root;
}

Code

daModesNode = daEditor("modes", "Soft prior: locked center mode (ring); click to add up to 5 more")
daMaskNode = daEditor("mask", "Hard mask: 8×8 blocks, starts upper-right 3×3; click blocks, then Save")
daModes = Generators.input(daModesNode)
daMask = Generators.input(daMaskNode)

Code

html`<div style="display:flex;gap:12px;flex-wrap:wrap;align-items:flex-start">
  <div style="flex:1 1 240px;max-width:260px">${daModesNode}</div>
  <div style="flex:1 1 240px;max-width:260px">${daMaskNode}</div>
</div>`

Code

{
  const {G, name, mirror, pts, ptCell, ghost} = daSecret;
  const N = daLadder[daN];
  const sigmaCells = daSigma, sig = sigmaCells / G;   // sd in unit-square coordinates
  const alpha0 = 0.5;                                  // shared base concentration
  const blue = "#2f64ad", amber = "#c2680c", slate = "#64748b";

  const centers = Array.from({length: G * G}, (_, i) => ({x: ((i % G) + .5) / G, y: (Math.floor(i / G) + .5) / G}));

  // per-cell densities -> per-cell Dirichlet alpha, normalized so mean(alpha)=alpha0
  const norm = arr => { const m = arr.reduce((s, v) => s + v, 0) / arr.length; return arr.map(v => alpha0 * v / (m || 1)); };
  const flatDen = centers.map(() => 1);
  const modeList = daModes.length ? daModes : [{c: G / 2, r: G / 2}];
  const mixDen = centers.map(ct => {
    let d = 0;
    for (const m of modeList) { const mx = (m.c + .5) / G, my = (m.r + .5) / G; d += Math.exp(-((ct.x - mx) ** 2 + (ct.y - my) ** 2) / (2 * sig * sig)); }
    return d / modeList.length + 1e-6;
  });
  const maskArr = daMask;                               // booleans, length G*G
  const alphaF = norm(flatDen);
  const alphaM = norm(mixDen);
  const KMask = maskArr.reduce((s, v) => s + (v ? 1 : 0), 0) || 1;
  const alphaK = maskArr.map(v => v ? alpha0 * (G * G) / KMask : 0);   // mean over allowed = alpha0

  // exact incremental Dirichlet-multinomial log-evidence on the shared point stream
  const ladder = daLadder;
  const evalModel = (alpha, allowNull) => {
    const A = alpha.reduce((s, v) => s + v, 0);
    const counts = new Float64Array(G * G);
    let lml = 0, nIn = 0, dead = 0;
    const rec = new Map([[0, 0]]);
    let li = 1;
    for (let n = 0; n < N; n++) {
      const c = ptCell[n];
      if (alpha[c] <= 0) { if (!dead) dead = n + 1; }
      else { lml += Math.log(alpha[c] + counts[c]) - Math.log(A + nIn); counts[c]++; nIn++; }
      if (ladder[li] === n + 1) { rec.set(n + 1, dead ? -Infinity : lml); li++; }
    }
    rec.set(N, dead ? -Infinity : lml);
    return {A, counts, dead, nIn, lml: dead ? -Infinity : lml,
            at: k => rec.has(k) ? rec.get(k) : NaN};
  };
  const mF = evalModel(alphaF), mM = evalModel(alphaM), mK = evalModel(alphaK);

  // 5000-point posterior-predictive cloud (Polya urn) at current N
  const cloud = (alpha, base, seed, M = 5000) => {
    const c2 = Float64Array.from(base.counts);
    const allowed = []; for (let i = 0; i < G * G; i++) if (alpha[i] > 0) allowed.push(i);
    let tot = base.A + base.nIn;
    const rng = mulberry32(seed);
    const out = [];
    for (let m = 0; m < M; m++) {
      let r = rng() * tot, cell = allowed[allowed.length - 1], acc = 0;
      for (const i of allowed) { acc += alpha[i] + c2[i]; if (acc >= r) { cell = i; break; } }
      c2[cell]++; tot++;
      out.push({x: ((cell % G) + rng()) / G, y: (Math.floor(cell / G) + rng()) / G});
    }
    return out;
  };

  // ---------- rendering ----------
  const box = (w, h, marks, title) => html`<div><div style="font-size:.78em;font-variant-caps:all-small-caps;letter-spacing:.04em;opacity:.65;margin:0 0 3px 4px">${title}</div>${themedPlot({width: w, height: h, x: {domain: [0, 1], axis: null}, y: {domain: [0, 1], axis: null}, marks: [Plot.frame({strokeOpacity: .2}), ...marks]})}</div>`;
  const W = 220, H = 220;

  const denCells = (alpha) => { const mx = Math.max(...alpha); return centers.map((ct, i) => ({x1: ct.x - .5 / G, x2: ct.x + .5 / G, y1: ct.y - .5 / G, y2: ct.y + .5 / G, v: alpha[i] / (mx || 1)})); };
  const denPanel = (alpha, scheme, title) => html`<div><div style="font-size:.78em;font-variant-caps:all-small-caps;letter-spacing:.04em;opacity:.65;margin:0 0 3px 4px">${title}</div>${themedPlot({width: W, height: H, x: {domain: [0, 1], axis: null}, y: {domain: [0, 1], axis: null}, color: {scheme, domain: [0, 1]}, marks: [Plot.rect(denCells(alpha), {x1: "x1", x2: "x2", y1: "y1", y2: "y2", fill: "v"})]})}</div>`;

  const reveal = daRevealTgl;
  const ghostLayer = reveal ? [Plot.rect(ghost, {x1: "x1", x2: "x2", y1: "y1", y2: "y2", fill: "currentColor", fillOpacity: .07})] : [];
  const obsShown = N <= 3000 ? pts.slice(0, N) : pts.slice(0, N).filter((_, i) => i % Math.ceil(N / 3000) === 0);
  const obs = box(W, H, [...ghostLayer, Plot.dot(obsShown, {x: "x", y: "y", r: 1.4, fill: "currentColor", fillOpacity: .5})],
    `Observed: ${N.toLocaleString()} point${N === 1 ? "" : "s"}${reveal ? ` - a ${mirror ? "mirrored " : ""}${name}` : ""}`);

  const race = html`<div style="grid-column: span 2"><div style="font-size:.78em;font-variant-caps:all-small-caps;letter-spacing:.04em;opacity:.65;margin:0 0 3px 4px">The evidence race - log evidence per observation</div>${(() => {
    const rows = [];
    for (const k of ladder) if (k > 0 && k <= Math.max(N, 1)) {
      const push = (m, nm) => { const v = m.at(k); if (isFinite(v)) rows.push({k: Math.log10(k), y: v / k, m: nm}); };
      push(mF, "flat"); push(mM, "soft mixture"); push(mK, "hard mask");
    }
    const dotRows = [];
    if (N > 0) { const push = (m, nm) => { const v = m.at(N); if (isFinite(v)) dotRows.push({k: Math.log10(N), y: v / N, m: nm}); }; push(mF, "flat"); push(mM, "soft mixture"); push(mK, "hard mask"); }
    const yv = rows.map(d => d.y); const lo = yv.length ? Math.min(...yv) : -3, hi = yv.length ? Math.max(...yv) : 0; const pad = .08 * (hi - lo + 1e-9);
    return themedPlot({width: 2 * W + 8, height: H, grid: true, x: {label: "N (log scale)", ticks: [0, 1, 2, 3, 4, 5], tickFormat: d => d >= 3 ? (10 ** d / 1000) + "k" : Math.round(10 ** d)}, y: {label: "log evidence / N", domain: [lo - pad, hi + pad], clamp: true},
      color: {domain: ["flat", "soft mixture", "hard mask"], range: [slate, blue, amber], legend: true},
      marks: [Plot.line(rows, {x: "k", y: "y", stroke: "m", strokeWidth: 2, curve: "monotone-x"}), Plot.dot(dotRows, {x: "k", y: "y", fill: "m", r: 5, stroke: "white", strokeWidth: 1.5})]});
  })()}</div>`;

  const postLbl = N === 0 ? "prior draws (no data yet)" : "5k draws";
  const postRow = [
    box(W, H, [...ghostLayer, Plot.dot(cloud(alphaF, mF, 111), {x: "x", y: "y", r: 1.1, fill: slate, fillOpacity: .5})], `flat ${postLbl}`),
    box(W, H, [...ghostLayer, Plot.dot(cloud(alphaM, mM, 222), {x: "x", y: "y", r: 1.1, fill: blue, fillOpacity: .5})], `soft ${postLbl}`),
    (mK.dead && N >= mK.dead)
      ? box(W, H, [...ghostLayer], "mask posterior - FALSIFIED")
      : box(W, H, [...ghostLayer, Plot.dot(cloud(alphaK, mK, 333), {x: "x", y: "y", r: 1.1, fill: amber, fillOpacity: .5})], `mask ${postLbl}`)
  ];

  const fmt = m => m.dead && N >= m.dead ? "−∞" : (m.at(N) === 0 && N === 0 ? "0.00" : (isFinite(m.at(N)) ? m.at(N).toFixed(1) : "—"));
  const readout = html`<p style="grid-column:1/-1;font-variant-numeric:tabular-nums;border-top:1px solid rgba(100,116,139,.2);padding-top:6px;font-size:.9em;margin-top:2px">
    ${N === 0 ? html`N = 0: every evidence is 0 (an empty dataset has probability 1 under any model).  The bottom row shows each prior generating from pure belief - flat static, your soft bump, your masked region.  Now add points and watch them learn.`
      : html`At N = ${N.toLocaleString()} - log evidence: flat <b style="color:${slate}">${fmt(mF)}</b>, soft mixture <b style="color:${blue}">${fmt(mM)}</b>, hard mask <b style="color:${amber}">${fmt(mK)}</b>${mK.dead && N >= mK.dead ? html` (a real point hit deleted cell at observation #${mK.dead.toLocaleString()})` : ""}.  Mask beats soft by <b>${(isFinite(mK.at(N)) && isFinite(mM.at(N))) ? (mK.at(N) - mM.at(N)).toFixed(1) : "—"}</b> nats - the price of its risk.`}
  </p>`;

  return html`<div style="display:grid;grid-template-columns:${W}px ${W}px ${W}px;gap:8px;align-items:start">
    ${denPanel(alphaF, "greys", "flat prior density")}
    ${denPanel(alphaM, "blues", "soft prior density")}
    ${denPanel(alphaK, "oranges", "mask prior density")}

    ${obs}
    ${race}

    ${postRow[0]}${postRow[1]}${postRow[2]}
    ${readout}
  </div>`;
}

Figure 14: A secret animal’s silhouette is the true density (reveal it any time). Three priors race on the SAME observed points: flat, a soft Gaussian mixture you place (a locked mode holds the center; click to add up to five more, width set by the σ slider), and a hard mask you carve in coarse blocks (starts as the upper-right 3×3; click blocks to add or delete, then Save). The density panels show each prior as the belief it is; the bottom panels show 5,000 draws from each posterior - the model’s picture of the animal. Drag N and watch: the flat prior crawls, the mask leaps ahead but risks -∞ the instant a real point lands in a deleted block, and the soft mixture - shaped with a little judgment - chases the mask down to the last few nats while never being able to die. That gap, small and safe versus large and risky, is Section 6’s whole trade. Drag N to 100,000 and watch even the flat prior finally resolve the animal - Lever 1, arriving late and expensive.

Read the three curves against each other, because they are Section 6’s postures made operational:

The flat prior cannot be confidently wrong: its only failure modes are blur and, at low concentration, parroting the observed points outright - Section 3’s mass collapse, drawn. It never generates the wrong animal, because it never believed anything specific enough to be that kind of wrong - and for exactly that reason it crawls, collecting its evidence one hard-won observation at a time. It reaches the truth eventually - drag \(N\) to 100,000 and watch it arrive - which is the promise and the price of Lever 1.
The mask is the opposite temperament. Carved well it leaps ahead - structure the data never had to buy - and holds its lead for as long as it lives. Carved with misplaced confidence it dies, the instant one real point lands in a cell you deleted, and its evidence is not large-and-negative but \(-\infty\) exactly, since the prior mass we placed on a point occurring in that location was already zero! Only a prior that can say impossible can be spectacularly wrong. But when the right choices are made (e.g. there are certain physical constraints that ensure that certain things happen), the truncated prior requires many fewer data points to learn the true underlying distribution.
And then the soft mixture, which is the whole reason this version exists. Watch where its curve sits: shaped with a little judgment - a few bumps roughly where a centered animal’s mass should be - it chases the mask down to the last handful of nats, and it never dies, because it kept a whisper of mass in every cell. The readout names that gap for you, and it is small. That is the honest shape of Section 6’s trade: the hard truncation’s advantage over a thoughtfully hedged prior is real but narrow, and it is bought with the risk of total collapse. The architect’s job is to decide, every time, whether those last few nats are worth the chance of being infinitely wrong - and the reason modern architectures truncate anyway is the subject of the next paragraph.

And notice the scale on which you just made that trade. Your mask was not one belief - it was \(64\) simultaneous per-cell declarations, most of them “impossible,” clicked in by hand, deleting most of the model’s prior volume in a minute of drawing. The impact here was somewhat small (but definitely noticeable). But in high dimensional settings like images, the impacts of these choices are enormous: a single click can delete billions of potential images from the model’s volume. Yes - ruling out something that could happen would lead your model to be spectacularly wrong if it did happen. But ruling out things that we know cannot happen using our human noggin’ is a must for ensuring that modern Generative AI models do not require billions upon billions of examples to “learn the state of the world”. Hard truncation is a high-stakes bet, but it is a bet that humans are very good at making.⁵

Two paths, one number, unequal bills. So there are two ways forward, and they raise the same objective. Collect more data and spend more compute; or do the hard human work of understanding the problem well enough to delete wasted volume. The mathematics prices them as substitutes. The world does not. One path’s costs are research labor, paid by the people doing the work. The other’s costs - electricity, water, land, the physical footprint of scale - are distributed to people who never joined the project. The choice between the levers is a real choice, made by someone, with consequences for someone else, and the field’s habit of treating it as methodologically neutral is itself a position.

The “wisp” has a name in the literature: the manifold hypothesis. See Bengio, Courville & Vincent, “Representation Learning” (2013); Fefferman, Mitter & Narayanan, “Testing the manifold hypothesis” (2016).

9.4 The catalogue of escapes

Here, then, is the taxonomy that organizes modern generative AI - and, not coincidentally, the rest of the content for this course. The hard object: the maximization we just stated, over an evidence we cannot evaluate. The field has found three escapes. Everything you have heard of is one of them.

Escape 1: Bound it. If you cannot compute the evidence, compute a floor under it, and raise the floor. That is the Evidence Lower BOund - the ELBO - and the machinery of variational inference, whose most famous child is the variational autoencoder. The cost of this escape is a gap: you optimize the floor, not the ceiling, and every compromise in the approximating family is a place the floor sags. (Where its compromises bite is exactly where the VQ-VAE written lecture begins.)

Escape 2: Decompose it in time. Break one impossible inference into a thousand easy ones: corrupt the data step by tiny step, and learn to undo each step - a chain of small conditional posteriors, each nearly Gaussian, each tractable. That is diffusion. The cost is a schedule: a thousand model calls to draw one sample, and a delicate choreography of noise levels.

Escape 3: Factor it exactly. For sequences, the chain rule is not an approximation. It is an identity:

\[ P(s_1, \dots, s_N) \;=\; \prod_{i=1}^{N} P(s_i \mid s_{<i}), \]

each factor a next-token predictive distribution. The cost is an ordering - you must decide what “previous” means - and a patience: one token at a time.

Now say out loud what Section 8 built. The probability pivot - predict the next thing given everything so far, score by log predictive likelihood, repeat - is autoregressive generative modeling. Train a network to output those predictive distributions, and the sequential validation scheme you watched grade ridge regression is the training objective of a GPT, verbatim: maximize \(\sum_i \log P(s_i \mid s_{<i})\). You did not follow a warm-up exercise two sections ago. You derived the loss function of the most consequential software artifact of the decade, using ridge regression and honesty.

One integral, three escapes. This figure is also the table of contents for the rest of the written lecture series.

Where each escape became an industry: Kingma & Welling, “Auto-Encoding Variational Bayes” (2013); Sohl-Dickstein et al. (2015) and Ho, Jain & Abbeel, “DDPM” (2020); Radford et al., “Language Models are Unsupervised Multitask Learners” (2019).

9.5 Coda: you were speaking Bayes the whole time

One last pass over familiar territory, with new eyes.

You called it	It is
weight decay	a Gaussian prior on weights
dropout	approximate averaging over an ensemble of sub-models
convolutional structure	a hard truncation of function space
data augmentation	an expansion of the prior over inputs
attention weights	a posterior over relevance
early stopping	implicit regularization of posterior geometry

Most of these are not analogies. They are identities - several exact, the rest principled approximations. Deep learning did not appear from nowhere; it emerged from decades of this reasoning, and it still speaks the language fluently even when its practitioners don’t notice they’re speaking it.

Srivastava et al., “Dropout” (JMLR 2014); Gal & Ghahramani, “Dropout as a Bayesian Approximation” (2016).

Where we stand - the final ledger. What we have: a language for reading every architecture as a commitment, every loss as a likelihood, every regularizer as a belief - and a compass, the marginal likelihood, with a precise meaning: honest one-step-ahead prediction, scored forever. What we lack: the machinery to compute any of it at the scale of images and language. Next: the machinery. First stop on this site: the VQ-VAE lecture - where Escape 1’s compromises meet a stupid simple fix.

Footnotes

The choice to refer to the marginal likelihood as encoding “surprise” is intentional since we used it to define the KL divergence above. The symmetry between the KL divergence of two distributions and the marginal likelihood of a model class is a deep one that will continue to appear throughout this series.↩︎
Note a particular notational stretch here using \(*\) as a superscript to denote the outcome at \(i\) as a random variable. The object returned here is a distribution!↩︎
Note that what we’re doing when we use \(K\)-fold CV is moving the uncertainty off of the parameter estimates themselves and instead inheriting uncertainty as a function of sampling variability. A reader with more statistical training may recognize this move as the same one that you learn about in probability classes when you discuss the difference between frequentist and Bayesian inference. Different day, same story with the same punchline that the two approaches are essentially equivalent in the limit!↩︎
This is actually a more nuanced claim than it seems on the surface. While there is general agreement that some metric of sequential likelihood form is the correct one for assessing generalization, which one is correct is still debated. The log marginal likelihood is the one grounded most closely in probability theory. The conditional log marginal likelihood (Lotfi et al. 2022) says that the LML doesn’t measure generalization as well as it could and proposes only beginning to compute the streaming predictive likelihood after some period of burn in. Other approaches like the WAIC (Watanabe 2010) are derived from expected properties in the singular model regime that is common in modern ML. All have their place. None answer every question. For supervised learning, any one of these could be viably used. For generative ML, the LML is the only one!↩︎
What does hard truncation look like in production generative AI? Convolutional locality and weight sharing are the canonical examples: a CNN declares that a pixel’s meaning cannot depend arbitrarily on far-away pixels, and that a feature detector working in one corner works everywhere - entire families of functions, deleted outright. Here is the subtle part: these truncations live in function space, not data space, so they cannot be falsified the way your mask can. No image receives zero probability from a CNN - feed it pure static and it computes a perfectly finite likelihood - so there is no \(-\infty\) moment, no loud death. If the locality assumption is wrong for a task, the failure is quiet: the true function is simply outside the model’s world, and the model fits the best wrong function available, forever, with no single observation able to object. Are the assumptions ever violated in practice? At the margins, maybe - vision transformers exist precisely to relax these truncations - and yet the generation-quality gap between a ViT and a well-built CNN is small while the CNN is far cheaper to train. That is the trade in the wild: a hard bet on locality that buys orders of magnitude in data and compute, priced against a failure mode that whispers instead of screams. And for completeness, production AI does contain truncations that die loudly, exactly like your mask - just not in the architecture. A language model’s fixed vocabulary assigns exactly zero probability to any string containing a symbol its tokenizer cannot represent: one out-of-vocabulary character, and the likelihood of the observed data is \(-\infty\), no appeal. Data-space truncations scream; function-space truncations whisper; the marginal likelihood is what hears both.↩︎

1. A model is a distribution

2. The fundamental problem

3. Under pressure: the collapse

4. Beliefs or degeneracy

5. The asymmetry: where priors actually matter

6. Real priors: truncate the hypothesis class

7. What the point estimate throws away

8. The probability pivot

9. Generative ML is Bayes with the training wheels off

9.1 The inner and outer problems

9.2 The integral is hard - and where cross-validation goes to die

9.3 What a generative model formally is

9.4 The catalogue of escapes

9.5 Coda: you were speaking Bayes the whole time

Further reading

Footnotes