Bayesian Statistics

A Bayesian model is a joint distribution over parameters and data, $p(\theta, x)$; everything else is conditioning. The course is a march of inference algorithms, deployed in order as exactness fails: conjugate exact → MAP + Laplace → MCMC (M–H, Gibbs, HMC) → variational. Each step trades exactness for reach, and the whole modern field lives at the last two stops — MCMC and VI are the load-bearing sections of this guide. The workflow wrapper is Box’s loop: build → infer → criticize → revise. Prerequisites cashed in: conjugacy and exp families (track 1 §2), posterior basics and MLE/MAP (track 3 §7), and the entire MCMC substrate — stationarity, detailed balance, ergodic theorem (track 5 §6).

1. The frame: posteriors and the proportionality workflow

• Bayes’ rule is the product and sum rules applied to the joint:

$p(\theta \mid x) = \frac{p(\theta, x)}{p(x)} = \frac{p(x \mid \theta)\,p(\theta)}{\int p(x \mid \theta)\,p(\theta)\,d\theta}$

Numerator: likelihood × prior — the model you wrote down. Denominator: the marginal likelihood / evidence — an integral over all of parameter space, and the villain of the entire course. Every algorithm below is a strategy for never computing it.

• The #1 working skill — recognize-the-shape: the denominator is constant in $\theta$, so

$p(\theta \mid x) \propto p(x \mid \theta)\,p(\theta)$

Multiply, simplify the $\theta$-dependence, match the shape to a known family; the normalizer comes free because that family is already normalized. You almost never integrate.

• The Gaussian instance of recognize-the-shape — complete the square. Any log-density that is quadratic in $z$ is a Gaussian, and you can read its parameters off the coefficients:

$\log p(z \mid \cdot) = -\tfrac{1}{2}z^TJz + h^Tz + \text{const} \;\Longrightarrow\; z \mid \cdot \sim N(J^{-1}h,\ J^{-1})$

$J$ (precision) and $h$ are the Gaussian’s natural parameters. Worked instance (factor analysis posterior — prior $N(0, I)$, likelihood $N(x \mid Wz, \Sigma)$): multiply, collect the quadratic and linear terms in $z$:

$-\tfrac{1}{2}z^Tz - \tfrac{1}{2}(x - Wz)^T\Sigma^{-1}(x - Wz) \;\Rightarrow\; J = I + W^T\Sigma^{-1}W, \qquad h = W^T\Sigma^{-1}x$

This one template computes every Gaussian conditional in the field — Gibbs conditionals, Kalman updates, GP posteriors.

• Improper priors (e.g. $p(\mu, \sigma^2) \propto \frac{1}{\sigma^2}$, the zero-pseudo-data limit of Normal-Inverse-Gamma): don’t integrate to 1, often still give a proper posterior — legal if you check. The cost arrives in §4: the evidence becomes undefined.
• Sequential = batch: today’s posterior is tomorrow’s prior; same answer as processing everything at once.
• Box’s loop: build the joint → infer → criticize (predictive checks, §9) → revise. Modeling is iteration, not specification.

Implications

• The log-joint $\log p(\theta, x)$ is the universal interface: MAP optimizers, M–H ratios, HMC dynamics, and the ELBO all consume only it (and its gradient). This is why probabilistic programming languages work — you write the joint, the inference engine does the rest.
• Recognize-the-shape is why conjugacy matters at all: it converts integration (hard, global) into pattern matching (easy, local).
• Complete-the-square is the reason “everything Gaussian is tractable” — linear-Gaussian models (regression, Kalman, factor analysis) stay closed under conditioning.

Core competency set

• Write Bayes’ rule from the joint; name what makes the denominator hard.
• Execute recognize-the-shape on a conjugate pair, and complete-the-square on a Gaussian product, reading off $J$ and $h$.
• Define improper priors, when they’re legal, and what they cost.

2. Conjugate exact inference and the hierarchical move

The exp-family machine from track 1 §2, driven in production:

• Likelihood of $n$ iid points: $\exp(\eta\sum t(x_i) - nA(\eta))$; conjugate prior $\propto \exp(\eta\phi - \nu A(\eta))$. Multiplying preserves the shape — updating is pseudo-data arithmetic:

$\phi \mathrel{\leftarrow} \phi + \sum_i t(x_i), \qquad \nu \mathrel{\leftarrow} \nu + n$

$\phi$ = pseudo-observations, $\nu$ = pseudo-count. Every conjugate-table row is this template.

• The Gaussian with both parameters unknown takes the Normal-Inverse-Gamma (≡ Normal-Inv-χ²) prior, $p(\mu, \sigma^2) = \mathrm{IGa}(\sigma^2 \mid \alpha, \beta)\,N(\mu \mid m, \sigma^2\lambda^{-1})$; the limit $\alpha = \beta = m = \lambda = 0$ recovers $\propto 1/\sigma^2$.
• Numbers (track 3’s formula, run once): $\sigma = 5$ known, prior $\mu \sim N(0, 10^2)$, $n = 4$, $\bar{x} = 8$. Precisions $\omega_0 = 0.01$, $n\omega = 0.16$:

$\mu \mid x \sim N\left(\frac{0.01(0) + 0.16(8)}{0.17},\ \frac{1}{0.17}\right) = N(7.53,\ 2.43^2)$

The prior, worth $\omega_0/\omega = 0.25$ observations, is outvoted 16-to-1.

The hierarchical move — eight schools, in full. Eight schools each report a treatment-effect estimate $\bar{y}_j$ with standard error $\sigma_j$. Three models:

1. No pooling: each school its own $\theta_j$, flat priors — estimates stay at $\bar{y}_j$, noisy schools stay noisy.
2. Complete pooling: one shared $\theta$ — as if one giant school; school identity erased.
3. Partial pooling (hierarchy): school effects drawn from a shared population:

$\bar{y}_j \mid \theta_j \sim N(\theta_j, \sigma_j^2), \qquad \theta_j \mid \mu, \tau \sim N(\mu, \tau^2), \qquad (\mu, \tau) \sim \text{hyperprior}$

The population level is itself learned from the data. Conditional on $(\mu, \tau)$, each school’s posterior mean is a precision-weighted compromise (complete-the-square again):

$E[\theta_j \mid \bar{y}_j, \mu, \tau] = \frac{\sigma_j^{-2}\,\bar{y}_j + \tau^{-2}\,\mu}{\sigma_j^{-2} + \tau^{-2}}$

• Read the limits: $\tau \to \infty$ ⟹ no pooling (the population prior is vague); $\tau \to 0$ ⟹ complete pooling (all schools forced to $\mu$). The data chooses $\tau$ — the data decides how much to pool.
• Numbers — the famous school A: $\bar{y}_A = 28$, $\sigma_A = 15$; suppose the population posterior centers at $\mu = 8$, $\tau = 6$. Precisions: $\sigma_A^{-2} = 0.0044$, $\tau^{-2} = 0.0278$:

$E[\theta_A \mid \cdot] = \frac{0.0044(28) + 0.0278(8)}{0.0322} \approx 10.8$

The headline effect of 28 shrinks to ~11 — most of it was noise, and the hierarchy knew because the other seven schools said effects that large don’t happen. Noisy groups borrow strength from the ensemble.

Implications

• Partial pooling is the Bayesian answer to regression-to-the-mean (track 4 §1) as a modeling device: shrinkage isn’t a correction you apply, it falls out of the joint.
• Hierarchies are conjugacy stacked: each level conditions on the one above; Gibbs sampling (§6) walks the levels.
• The same structure powers multilevel regression, random effects (track 4’s mixed models), and empirical Bayes; “how much do groups share?” becomes a parameter ($\tau$) instead of a modeling fiat.

Core competency set

• Run the pseudo-data template generically and for Normal-Inverse-Gamma; compute the Normal–Normal numeric.
• Write the eight-schools hierarchy, derive the precision-weighted shrinkage, read both $\tau$ limits, and run the school-A numbers.
• Explain “borrowing strength” in one sentence.

3. Bayesian linear regression — and ridge unmasked

• Model: $y_n \mid x_n, w, \sigma^2 \sim N(w^Tx_n, \sigma^2)$; conjugate prior $\text{Inv-}\chi^2(\sigma^2 \mid \nu, \tau^2)\,N(w \mid \mu, \sigma^2\Lambda^{-1})$.
• Posterior by recognize-the-shape — same family, sufficient statistics absorbed:

$\Lambda' = \Lambda + \sum_n x_nx_n^T, \qquad \mu' = \Lambda'^{-1}\Big(\Lambda\mu + \sum_n y_nx_n\Big), \qquad \nu' = \nu + N$

Pseudo-data again: the prior contributes $\Lambda$ fake $xx^T$‘s and $\Lambda\mu$ fake $yx$‘s.

• Uninformative limit ($\nu, \Lambda \to 0$): MAP = MLE = OLS — track 4 rung 2, met from the other side.
• Ridge unmasked. Set $\Lambda = \lambda I$, $\mu = 0$:

$w_{MAP} = \Big(\lambda I + \sum_n x_nx_n^T\Big)^{-1}\sum_n y_nx_n = (X^TX + \lambda I)^{-1}X^Ty$

Ridge is the posterior mode under a spherical Gaussian prior — the penalty is the prior, $\lambda$ is a prior precision. Track 4 §8’s “targeted variance repair” and “shrink toward the prior mean” are one operation in two vocabularies. Lasso = MAP under a Laplace (double-exponential) prior — its sharp peak at zero is what produces exact zeros.

• Posterior predictive — integrate the parameters out, never plug in:

$p(y_{N+1} \mid x_{N+1}, \text{data}) = \int p(y_{N+1} \mid x_{N+1}, w, \sigma^2)\,p(w, \sigma^2 \mid \text{data})\,dw\,d\sigma^2$

Here analytic: Gaussian × Inv-χ² ⟹ Student-t — track 1 §8’s “normal with a noisy scale,” rediscovered as parameter uncertainty. Credible bands flare outside the data; in general estimate by Monte Carlo (sample parameters, average conditional densities).

Implications

• Every regularizer is a prior. The penalty-as-prior dictionary (ridge↔Gaussian, lasso↔Laplace) means regularization tuning is prior elicitation — and exposes what assumptions a penalty silently makes.
• Prediction with integrated-out parameters is systematically better-calibrated than plug-in: the t’s heavy tails are honesty about not knowing $\sigma$.
• Frequentist answers reappear as limits (flat priors, modes) — the frameworks disagree at small $n$ and informative priors, exactly where the prior has something to say.

Core competency set

• Derive the posterior updates and narrate them as pseudo-data.
• Reproduce ridge-as-MAP in two lines; state the lasso analogue.
• Explain why the predictive is a t and what integrating out buys over plugging in.

4. The marginal likelihood: evidence and Occam

The villain rehabilitated as a model-comparison tool:

• $p(x) = \int p(x \mid \theta)\,p(\theta)\,d\theta = E_{p(\theta)}[\,p(x \mid \theta)\,]$ — the probability the prior-averaged model assigns the data. Not best-case fit; expected fit, before seeing the data.
• The Z-ratio trick, displayed once (works in any conjugate family). Write the unnormalized posterior and multiply-divide by its known normalizer $Z'$:

$p(x) = \int p(x \mid \theta)\,p(\theta)\,d\theta = \frac{c}{Z}\int \tilde{p}(\theta \mid x)\,d\theta = \frac{c}{Z}\cdot Z'\underbrace{\int \frac{\tilde{p}(\theta \mid x)}{Z'}\,d\theta}_{=\,1} = c\,\frac{Z'}{Z}$

($c$ collects likelihood constants like $(2\pi)^{-N/2}$.) Evidence = ratio of posterior to prior normalizing constants — the same multiply-by-one move that proves conjugate predictives.

• Automatic Occam’s razor. In polynomial-degree selection the evidence factors into a fight:
  • fit: $\tau'^2 \propto y^T(I - H)y$ (residual error) falls with degree ⟹ pushes evidence up;
  • complexity: $\frac{|\Lambda|^{1/2}}{|\Lambda'|^{1/2}}$ falls as dimension grows — more parameters spread the prior’s mass thinner over data space, so any particular dataset receives less of it. The maximum sits at “complex enough.” No penalty was added by hand — it lives in the integral.
• Caveats: undefined for improper priors (infinitely spread prior gives every dataset measure zero), and genuinely sensitive to prior choices.

Implications

• Ratios of evidences are Bayes factors — the Bayesian replacement for hypothesis tests; Occam is built in rather than bolted on (contrast $C_p$/AIC’s additive corrections, track 4 §7).
• Laplace-approximating this integral (§5) and keeping the $O(\log n)$ terms yields BIC — the cheap shadow of the evidence.
• The prior-sensitivity caveat is what pushes practice toward the predictive currency (§9) when priors are diffuse.

Core competency set

• Define the evidence as a prior expectation; run the Z-ratio manipulation.
• Tell the Occam story with both factors named.
• State the caveats and what each pushes you toward.

5. GLMs, MAP, and the Laplace approximation

First failure of conjugacy — non-Gaussian observations:

• GLM recipe: observations exp-family, $p(y_n \mid \eta_n) = h(y_n)\exp(\langle t(y_n), \eta_n\rangle - A(\eta_n))$, with $\eta_n = w^Tx_n$. The mean is $E[t(y)] = \nabla A(\eta)$ (track 1’s identity), so the canonical mean function is $\nabla A$ — logistic for Bernoulli, $\exp$ for Poisson. The link isn’t a modeling choice; it’s the family’s own gradient.
• Logistic regression + Gaussian prior $N(0, \sigma^2I)$: the posterior has no closed form (the $1 + e^{w^Tx_n}$ denominators never collapse). Retreat to the mode:

$\nabla_w\mathscr{L}(w) = (y - \hat{y})^TX - \frac{1}{\sigma^2}w, \qquad \hat{y} = f(Xw)$

— an error signal against the covariates, minus the prior’s pull to zero. The Hessian is negative definite ($f' \geq 0$, plus the prior), so the log-posterior is concave: Newton/IRLS finds the unique MAP.

• Laplace approximation — promote the mode to a distribution. Second-order Taylor at $w_{MAP}$ (linear term dead there):

$\log p(w \mid x) \approx \log p(w_{MAP} \mid x) - \tfrac{1}{2}(w - w_{MAP})^T\big[-\nabla^2\mathscr{L}\big](w - w_{MAP})$

$\Rightarrow\quad p(w \mid x) \approx N\Big(w_{MAP},\ \big[-\nabla^2\mathscr{L}(w_{MAP})\big]^{-1}\Big)$

Picture to hold: a skewed posterior bump with a symmetric Gaussian snugged to its summit — exact at the top, wrong in the tails.

• Bernstein–von Mises: as $n \to \infty$ the posterior is asymptotically normal, centered at the MLE with variance $[-\ell'']^{-1} \approx \frac{1}{nI}$ — track 3 §7’s large-sample posterior with its proper name; the log-likelihood Hessian converges to Fisher information, one observation’s curvature at a time.

Implications

• Laplace gives cheap uncertainty for anything you can MAP-train — including neural networks (Laplace on the last layer is a modern serving trick). Cost: one Hessian.
• BvM is the peace treaty: with enough data, Bayesian and frequentist uncertainty coincide. All the disagreement lives at small $n$ — which is where you actually need the prior.
• Concavity of exp-family GLM posteriors is why logistic/Poisson regression are reliable workhorses: one optimum, no restarts.

Core competency set

• Build a GLM and explain why the canonical link is $\nabla A$.
• Write the logistic MAP gradient; argue concavity.
• Derive Laplace from the Taylor expansion; state BvM and its track-3 identity.

6. MCMC — sampling the posterior

When no approximation is trusted, sample. Track 5 §6 built the theory; here it earns its living. This and §8 are the two sections to over-learn.

Why Monte Carlo at all. Everything you want is an expectation — posterior means, event probabilities, predictives are all $E_{p(\theta \mid x)}[f(\theta)]$. Quadrature’s cost grows exponentially in dimension; the MC estimator $\hat{f} = \frac{1}{S}\sum_s f(\theta^{(s)})$ is unbiased with standard error $O(S^{-1/2})$ independent of dimension. High dimensions are exactly where Bayes lives, so MC wins there.

The obstruction and the escape. Sampling $p(\theta \mid x)$ directly is as hard as computing $p(x)$. But ratios of posteriors need only the joint:

$\frac{p(\theta \mid x)}{p(\theta' \mid x)} = \frac{p(\theta, x)/p(x)}{p(\theta', x)/p(x)} = \frac{p(\theta, x)}{p(\theta', x)}$

So build a Markov chain that only ever consults ratios — then the evidence never appears.

Metropolis–Hastings, derived from detailed balance. We want a transition kernel $\pi(\theta \mid \theta')$ with the posterior as its stationary distribution; detailed balance (track 5 §5) suffices. Write the kernel as propose-then-accept, $\pi(\theta \mid \theta') = q(\theta \mid \theta')\,a(\theta' \to \theta)$, and impose DBC:

$p(\theta' \mid x)\,q(\theta \mid \theta')\,a(\theta' \to \theta) = p(\theta \mid x)\,q(\theta' \mid \theta)\,a(\theta \to \theta')$

$\Rightarrow\quad \frac{a(\theta' \to \theta)}{a(\theta \to \theta')} = \frac{p(\theta, x)\,q(\theta' \mid \theta)}{p(\theta', x)\,q(\theta \mid \theta')}$

Only the ratio of acceptance probabilities is pinned down — so make them as large as possible by setting the bigger one to 1:

$a(\theta' \to \theta) = \min\left\{1,\ \frac{p(\theta, x)\,q(\theta' \mid \theta)}{p(\theta', x)\,q(\theta \mid \theta')}\right\}$

• Symmetric $q$ ⟹ accept uphill always, downhill with probability = the density ratio: stochastic hill climbing that provably samples the right thing. Any $q$ that can reach everywhere gives an ergodic chain; the ergodic theorem (track 5) then makes time averages converge to posterior expectations.
• Designing $q$ is the whole game: too-small steps ⟹ high acceptance but snail mixing; too-big ⟹ constant rejection. (Optimal random-walk acceptance ≈ 0.23 in high dimension — you want rejections.)

Gibbs sampling = M–H whose proposal is the full conditional, $q(\theta \mid \theta') = p(\theta_j \mid \theta'_{-j}, x)$, changing one block at a time. The acceptance ratio cancels identically — display it once, with $\theta_{-j} = \theta'_{-j}$ fixed:

$\frac{p(\theta, x)\,q(\theta' \mid \theta)}{p(\theta', x)\,q(\theta \mid \theta')} = \frac{p(\theta_j \mid \theta_{-j}, x)\,p(\theta_{-j}, x)\cdot p(\theta'_j \mid \theta_{-j}, x)}{p(\theta'_j \mid \theta_{-j}, x)\,p(\theta_{-j}, x)\cdot p(\theta_j \mid \theta_{-j}, x)} = 1$

• Acceptance probability 1, always. Cycle through coordinates/blocks for ergodicity. Gibbs shines in conditionally conjugate models — each conditional is a known family (recognize-the-shape / complete-the-square supply them).
• Where Gibbs dies: correlated coordinates. The picture to hold: axis-parallel zigzag steps trapped inside a long tilted ellipse — each conditional slice is tiny, so the chain inches along the ridge. Note it’s the correlation, not the elongation: an axis-aligned ellipse samples perfectly.
• Fixes: reparametrize to decorrelate; block correlated coordinates and sample them jointly; collapse — integrate a parameter out analytically (exp-family conjugacy makes this possible) and run the chain on the lower-dimensional marginal. Collapsing trades parallelism for mixing.

Diagnostics — the estimator is consistent but correlated and transiently biased.

• Trace plots: flat stretches = rejections; parameters mix at different rates. Burn-in: discard early samples still under the initialization’s influence (bias decays faster than variance).
• Autocorrelation discounts your sample size. Effective sample size:

$S_{\mathrm{eff}} = \frac{S}{1 + 2\sum_{\ell=1}^\infty \mathrm{acf}(\ell)}$

• Benchmark it on AR(1), $\theta_s = \rho\theta_{s-1} + \epsilon_s$, where the whole computation is a geometric series: $\mathrm{acf}(\ell) = \rho^\ell$, so

$1 + 2\sum_{\ell \geq 1}\rho^\ell = 1 + \frac{2\rho}{1-\rho} = \frac{1+\rho}{1-\rho} \quad\Rightarrow\quad S_{\mathrm{eff}} = S\,\frac{1-\rho}{1+\rho}$

Numbers: $\rho = 0.9$ ⟹ $S_{\mathrm{eff}} = S/19$ — a 100k chain is worth ~5k independent draws. This number is why the gradient methods below exist.

Use the gradient. Autodiff makes $\nabla_\theta \log p(\theta, x)$ essentially free, and it changes the game:

• MALA: drift the proposal uphill, $q(\theta' \mid \theta) = N\big(\theta + \tau\nabla\log p(\theta, x),\ 2\tau^2I\big)$. Asymmetric, so both directions’ $q$‘s enter the acceptance ratio. Cuts the steps-per-independent-sample from $O(d^2)$ toward $O(d^{4/3})$.
• HMC — the big idea: stop diffusing, start coasting. Augment with a momentum variable $p$ and define a physical system:

$H(q, p) = \underbrace{U(q)}_{-\log p(\theta = q,\ x)} + \underbrace{K(p)}_{\frac{1}{2}p^TM^{-1}p}, \qquad p(q, p) \propto e^{-H(q,p)} = p(\theta \mid x) \times N(p \mid 0, M)$

The joint factorizes — sample momenta fresh, simulate the dynamics, keep only $q$, and the $q$-marginal is exactly the posterior. The dynamics follow Hamilton’s equations ($\dot{q} = M^{-1}p$, $\dot{p} = -\nabla U$): a frictionless puck coasting over the negative log posterior landscape — it converts potential to kinetic energy and back, traveling far along level sets instead of dithering. Energy is conserved ⟹ proposals from exact dynamics would be accepted with probability 1.

• Discretization must respect the geometry. The leapfrog integrator alternates half-steps:

$p_{t+\epsilon/2} = p_t - \tfrac{\epsilon}{2}\nabla U(q_t), \qquad q_{t+\epsilon} = q_t + \epsilon\,M^{-1}p_{t+\epsilon/2}, \qquad p_{t+\epsilon} = p_{t+\epsilon/2} - \tfrac{\epsilon}{2}\nabla U(q_{t+\epsilon})$

Each update is a **shear** (one variable moves, by an amount depending only on the *other*) — determinant 1, volume preserved — which is exactly what forward Euler lacks (its trajectories spiral out). The small leftover energy error is repaired by one M–H accept/reject.

• NUTS auto-tunes the trajectory length: simulate until the path starts doubling back on itself (the U-turn), avoiding both wasted laps and premature stops. NUTS-on-HMC is the default engine of Stan and PyMC — when someone “fits a Bayesian model” today, this is usually what runs.
• The practical art: mix and match — Gibbs the conditionally conjugate coordinates, HMC the rest; collapse what you can.

Implications

• MCMC’s contract: asymptotically exact — the gold standard against which approximations (§8) are judged; its price is correlation (ESS) and per-step cost.
• The gradient era reorganized the field: HMC/NUTS made “write the joint, autodiff it, sample” a button — modeling effort shifted from designing samplers to designing models (Box’s loop sped up).
• Diagnostics are not optional bureaucracy: an unconverged chain looks like a converged one locally. Trace + ESS (+ multi-chain $\hat{R}$ in practice) are the seatbelt.

Core competency set

• Explain dimension-free MC error and the ratio escape from the evidence.
• Derive M–H acceptance from detailed balance; display the Gibbs cancellation; draw the correlated-ellipse failure and name the three fixes.
• Run the AR(1) ESS derivation and the $\rho = 0.9$ number.
• Tell the HMC story: augmentation, factorized joint, coasting, energy conservation, leapfrog-as-shears, NUTS.

7. Latent variable models: mixtures and mixed membership

• Bayesian mixture model — the generic clustering joint: $\pi \sim \mathrm{Dir}(\alpha)$; components $\eta_k \sim$ conjugate prior; assignments $z_n \sim \mathrm{Cat}(\pi)$; data $x_n \mid z_n \sim \text{expfam}(\eta_{z_n})$. The joint factors — exactly what Gibbs wants.
• The Gibbs sweep, every conditional in closed form:

$z_n \mid \cdot \;\sim\; \mathrm{Cat}\Big(\tfrac{\pi_k\,p(x_n \mid \eta_k)}{\sum_j \pi_j\,p(x_n \mid \eta_j)}\Big) \quad\text{(the responsibilities)}$

$\eta_k \mid \cdot \;\propto\; \exp\Big(\eta_k\big(\phi + \textstyle\sum_{n: z_n = k}t(x_n)\big) - (\nu + N_k)A(\eta_k)\Big) \quad\text{(conjugate update, component-}k\text{ data)}$

$\pi \mid \cdot \;\sim\; \mathrm{Dir}(\alpha_1 + N_1, \ldots, \alpha_K + N_K) \quad\text{(pseudo-counts + cluster counts)}$

• Degenerate limits locate the classics: MAP coordinate ascent on $z$ with shared spherical Gaussians ⟹ k-means; $N_k \to \infty$ collapses conditionals to point estimates ⟹ EM-like behavior. Soft assignment is the general case; hard assignment is its zero-temperature limit.
• Mixed membership / LDA: grouped data — each document mixes the same $K$ topics with its own proportions. Sharing at the corpus level, variability at the document level. The posterior’s tension — explain each document with few topics vs each topic with few words — is what makes topics sharp. Exchangeability (de Finetti: any exchangeable joint is a mixture of iid models) is the license for bag-of-words.

Implications

• The mixture Gibbs sweep is the template for all conditionally conjugate latent-variable models — HMMs (transitions linked $z$‘s), topic models, factor models all reuse it with one conditional swapped.
• Responsibilities are the bridge to EM (§8): E-step = expected responsibilities, M-step = their weighted conjugate updates.

Core competency set

• Write the mixture joint and all three Gibbs conditionals from memory.
• Place k-means as a degenerate limit; narrate the LDA tension and the de Finetti license.

8. Variational inference — optimization replaces sampling

MCMC is asymptotically exact but slow to decorrelate; VI trades bias for variance: pick a tractable family $q(\theta; \lambda)$, then optimize $\lambda$ to make $q$ close to the posterior. You give up exactness-in-the-limit and get back speed, determinism, and SGD-scale data handling. Modern ML’s generative stack lives here.

• Closeness = KL divergence, and the direction is a real decision:

$\mathrm{KL}(q \| p) = E_q\Big[\log\frac{q(\theta)}{p(\theta \mid x)}\Big]$

Nonnegative, zero iff equal, an expectation under $q$ — which we chose to be easy. The direction makes $q$ mode-seeking / zero-forcing: where $p \approx 0$, any $q$-mass costs infinitely, so $q$ retreats inside one mode and underestimates spread. (The reverse KL, $E_p[\log p/q]$, is mass-covering — but you can’t take expectations under the posterior you don’t have.) Picture to hold: a bimodal posterior with $q$ locked snugly onto one mode, ignoring the other entirely.

• The ELBO derivation — removing the unknown posterior from the objective. Expand $p(\theta \mid x) = p(\theta, x)/p(x)$ inside the KL:

$\mathrm{KL}(q \| p) = E_q[\log q(\theta; \lambda)] - E_q[\log p(\theta, x)] + \log p(x)$

$\Rightarrow\quad \log p(x) = \underbrace{E_q[\log p(\theta, x)] - E_q[\log q(\theta; \lambda)]}_{\mathrm{ELBO}\ \mathscr{L}(\lambda)} + \underbrace{\mathrm{KL}(q \| p)}_{\geq 0}$

The evidence is constant in $\lambda$: maximizing the ELBO ⟺ minimizing the KL, and the ELBO needs only the joint (the universal interface again). Picture to hold: a fixed bar of height $\log p(x)$ split into ELBO below, KL above — pushing the ELBO up squeezes the KL out. And since $\mathrm{KL} \geq 0$, the ELBO really is a lower bound on the evidence — usable for model comparison too.

• Two readings of the same objective, each guarding against a failure:

$\mathscr{L}(\lambda) = E_q[\log p(x \mid \theta)] - \mathrm{KL}\big(q \,\|\, p(\theta)\big) \qquad \text{(fit the data, stay near the prior)}$

$\mathscr{L}(\lambda) = E_q[\log p(\theta, x)] + H[q] \qquad \text{(seek the joint's mass, but keep entropy — don't collapse to the MAP point)}$

• Mean field: $q(\theta; \lambda) = \prod_j q(\theta_j; \lambda_j)$ — all posterior dependence discarded by fiat. What’s lost: correlations, and (via zero-forcing) variances come out too small. What’s gained: every expectation factorizes.
• CAVI — coordinate ascent on the ELBO. Hold $q_{-j}$ fixed; as a function of $q_j$ the ELBO is $-\mathrm{KL}(q_j \| \tilde{p}_j) + \text{const}$, maximized by

$q^*_j(\theta_j) \;\propto\; \exp\Big\{E_{q_{-j}}\big[\log p(\theta_j \mid \theta_{-j}, x)\big]\Big\}$

CAVI is Gibbs with expectations in place of samples: same conditional structure, same conjugate algebra, but a deterministic sweep that monotonically increases the ELBO. (For LDA: same updates as the Gibbs sampler, expectations swapped in.)

• EM placed exactly. With point-estimated parameters $\Theta$ and latent variables $z$: the ELBO satisfies

$\mathscr{L}(q, \Theta) = \log p(x \mid \Theta) - \mathrm{KL}\big(q(z)\,\|\,p(z \mid x, \Theta)\big)$

so the bound is tight exactly when $q$ = the true posterior over latents — that’s the E-step; the M-step maximizes over $\Theta$. EM = alternating exact-fit of $q$ and parameter maximization. Variational EM = EM where the E-step posterior is itself out of reach and must be approximated within a family — opening the approximation gap (best-in-family KL > 0).

• Stochastic / black-box VI — what makes VI scale. Optimize the ELBO by SGD; the obstacle is $\nabla_\lambda E_{q(z;\lambda)}[f(z)]$ (the parameter sits inside the sampling distribution). The reparameterization trick moves it:

$z = g(\lambda, \epsilon) = \mu + \sigma\epsilon,\ \ \epsilon \sim N(0, I) \quad\Rightarrow\quad \nabla_\lambda E_{q(z;\lambda)}[f(z)] = E_\epsilon\big[\nabla_\lambda f(g(\lambda, \epsilon))\big]$

Sample $\epsilon$‘s, differentiate through $g$ — a Monte Carlo gradient with low variance, requiring $f$, the density, and $g$ differentiable. (Discrete latents need the score-function/REINFORCE estimator instead — unbiased but high-variance.)

• VAEs — amortize. Instead of optimizing a separate $\lambda_n$ per data point, learn an inference network $x_n \mapsto \lambda_n$ (the encoder) jointly with the generative model (the decoder), training both by reparameterized SGD on the ELBO. Two gaps now stack: the approximation gap (family can’t reach the posterior) + the amortization gap (the network doesn’t even find the family’s best member). Sanity anchor: for linear factor analysis the optimal encoder is exactly the analytic posterior map $J^{-1}h$ from §1 — a pseudoinverse-like linear function; the network’s job is to learn the decoder’s inverse.

Implications

• The modern division of labor: MCMC when you need fidelity (science, small-to-medium models, calibrated uncertainty), VI when you need scale or speed (embedded in larger systems, millions of data points, deep generative models). The bias–variance trade is the organizing axis of contemporary inference.
• Standing warning: mean-field VI underestimates posterior variance — never read calibrated uncertainty off a vanilla mean-field fit.
• The ELBO is the loss function of deep generative modeling: VAEs and their descendants are this section with neural networks for $f$ — diffusion models train on a (hierarchical) ELBO too.
• CAVI ⟷ Gibbs is a two-way street: a model you can Gibbs you can CAVI, and vice versa — derive one, get the other’s structure free.

Core competency set

• Derive the ELBO from the KL; state evidence = ELBO + KL with the bar picture; explain why ELBO-maximization needs only the joint.
• Explain mode-seeking from the KL direction and its underdispersion consequence, with the bimodal picture.
• Write the CAVI update and its Gibbs correspondence; place EM, variational EM, and VAEs with their two gaps.
• Derive the reparameterization trick and say what breaks without it.

9. Model criticism — closing Box’s loop

• Posterior predictive checks: simulate replicate datasets from the posterior predictive; compare a test statistic $T$ between replicates and observed data; the tail fraction is a posterior predictive p-value. Internal coherence: “does my fitted model generate data that looks like mine?”
• Predictive scoring (the external currency): log predictive density on held-out data, $\log\int p(y' \mid \theta)\,p(\theta \mid y)\,d\theta$ — integrate over the posterior, never plug in. CV is the standard; WAIC/DIC approximate it from the posterior itself (AIC matches only when MAP ≈ MLE).
• Sensitivity analysis: perturb priors and structure, watch the conclusions; extrapolations are fragile, interpolations robust.
• The two currencies divide cleanly: evidence compares models as priors (§4); prediction compares them as forecasters. Diffuse or improper priors ⟹ only the second is available.

Core competency set

• Run a PPC in words: statistic, replicates, tail probability.
• Contrast evidence vs predictive scoring and when each is the right currency.

10. Memorize cold

The instant-recall layer — these should cost nothing to produce; the fluency flashcards drill exactly this list.

• Bayes from the joint; posterior ∝ likelihood × prior; recognize-the-shape workflow.
• Complete-the-square: $\log p(z) = -\frac{1}{2}z^TJz + h^Tz + c \Rightarrow z \sim N(J^{-1}h, J^{-1})$; factor-analysis instance $J = I + W^T\Sigma^{-1}W$, $h = W^T\Sigma^{-1}x$.
• Conjugate template: $\phi \mathrel{\leftarrow} \phi + \sum t(x_i)$, $\nu \mathrel{\leftarrow} \nu + n$; Normal-Inverse-Gamma for $(\mu, \sigma^2)$; improper limit $\propto 1/\sigma^2$.
• Eight schools: $\theta_j \sim N(\mu, \tau^2)$; shrinkage $E[\theta_j \mid \cdot] = \frac{\sigma_j^{-2}\bar{y}_j + \tau^{-2}\mu}{\sigma_j^{-2} + \tau^{-2}}$; $\tau \to 0$ complete pooling, $\tau \to \infty$ none.
• Regression updates $\Lambda' = \Lambda + \sum x_nx_n^T$, $\mu' = \Lambda'^{-1}(\Lambda\mu + \sum y_nx_n)$; ridge = MAP under Gaussian prior, lasso = MAP under Laplace prior.
• Predictive: integrate parameters out; Gaussian × Inv-χ² ⟹ Student-t.
• Evidence $= E_{prior}[\text{lik}] = c\,Z'/Z$; Occam = fit ↑ vs prior-mass-spreading ↓; Bayes factors; undefined for improper priors; BIC = Laplace’s skeleton.
• Canonical link = $\nabla A$; logistic MAP gradient $(y - \hat{y})^TX - w/\sigma^2$; concave.
• Laplace: $N(\theta_{MAP}, [-\nabla^2\log p]^{-1})$; Bernstein–von Mises → $N(\hat\theta_{MLE}, \frac{1}{nI})$.
• M–H from DBC: $a = \min\{1, \frac{p(\theta,x)q(\theta'|\theta)}{p(\theta',x)q(\theta|\theta')}\}$; ratios kill the evidence; ~0.23 optimal RW acceptance.
• Gibbs = full-conditional proposal ⟹ acceptance 1 (the cancellation); fails on correlated coordinates; fixes: reparametrize / block / collapse.
• ESS $= S/(1 + 2\sum\mathrm{acf})$; AR(1): $\mathrm{acf}(\ell) = \rho^\ell \Rightarrow S\frac{1-\rho}{1+\rho}$; $\rho = 0.9 \Rightarrow S/19$.
• HMC: $H = U + K$, $U = -\log p(\theta, x)$; joint factorizes so the $q$-marginal is the posterior; leapfrog = volume-preserving shears; energy conservation ⟹ near-1 acceptance; NUTS stops at the U-turn.
• Mixture Gibbs: responsibilities / conjugate $\eta_k$ / $\mathrm{Dir}(\alpha + N_k)$; k-means = degenerate limit.
• $\log p(x) = \mathrm{ELBO} + \mathrm{KL}(q \| p)$; ELBO $= E_q[\log p(\theta, x)] + H[q] = E_q[\log p(x \mid \theta)] - \mathrm{KL}(q \| \text{prior})$.
• KL$(q\|p)$ is mode-seeking/zero-forcing ⟹ mean-field VI underestimates variance.
• CAVI: $q_j^* \propto \exp E_{-j}[\log p(\theta_j \mid \theta_{-j}, x)]$ — Gibbs with expectations.
• EM: bound tight ⟺ $q$ = latent posterior; variational EM adds the approximation gap; VAEs add the amortization gap.
• Reparameterization: $z = \mu + \sigma\epsilon \Rightarrow \nabla_\lambda E_q[f] = E_\epsilon[\nabla_\lambda f]$.
• Named moves: recognize-the-shape; complete-the-square; pseudo-data update; integrate-out; multiply-by-$Z'/Z$; ratios-kill-the-evidence; trade-bias-for-variance; augment-then-marginalize (HMC’s momentum).