Skip to content
Four Roads to a Posterior, home

MAST90125 · Bayesian Statistical Learning · 2023

Four roads to a posterior

One small dose-finding study, one Bayesian logistic regression, and four ways to compute its posterior: a Laplace approximation, Metropolis-Hastings, Hamiltonian Monte Carlo and expectation propagation. The original R code now runs in your browser and reproduces the 2023 output draw for draw.

Open the explorer Watch HMC moveStep through EP
New here? Take the guided tour: three short walkthroughs with transcripts
Posterior of the dose slope β₁corrected · seed 90125
0.00.20.40.6no dose effect
  • Laplace
  • MH
  • HMC
  • EP

The question

Does a bigger dose help more patients?

A trial of an experimental influenza medication gave one of seven doses (0 to 6) to ten patients each, and a week later recorded whether their symptoms had improved. The assignment asked for a logistic regression of improvement on dose with a flat prior, fitted four different ways, and a comparison of what each method says about the coefficients.

It is a deliberately tiny problem: seventy patients and two parameters. That makes it ideal for seeing how the methods differ, because the true posterior can be computed on a grid and drawn underneath all of them.

Patients improved out of ten at each dose
dose0123456
improved3544768
yi∼Bin⁡(10, pi)logit⁡pi=β0+β1xip(β)∝1\begin{gathered} y_i \sim \operatorname{Bin}(10,\, p_i) \\ \operatorname{logit} p_i = \beta_0 + \beta_1 x_i \\ p(\beta) \propto 1 \end{gathered}

37 of 70 patients improved overall, rising from 3 of 10 at dose 0 to 8 of 10 at dose 6.

The four roads

Two approximations, two samplers

Each method is a line-by-line port of the function in the submitted R Markdown, which in turn adapted the MAST90125 lecture code.

  1. Laplace approximation

    01

    approximation

    p(β∣y)≈N ⁣(β^, (X⊤WX)−1)p(\beta \mid y) \approx N\!\big(\hat\beta,\ (X^\top W X)^{-1}\big)

    Find the posterior mode and fit a Gaussian to the curvature there. With a flat prior this is exactly R's glm(), fitted by iteratively reweighted least squares.

  2. Metropolis-Hastings

    02

    sampler

    β∗∼N ⁣(β(t), c2Σ^),  α=min⁡ ⁣(1,p(β∗∣y)p(β(t)∣y))\beta^\ast \sim N\!\big(\beta^{(t)},\ c^2 \hat\Sigma\big),\ \ \alpha = \min\!\Big(1, \tfrac{p(\beta^\ast \mid y)}{p(\beta^{(t)} \mid y)}\Big)

    Propose a random jump scaled from the GLM covariance and accept it with the posterior ratio. Simple and exact in the limit, but successive draws are strongly correlated.

  3. Hamiltonian Monte Carlo

    03

    sampler

    H(β,ϕ)=−log⁡p(β∣y)+12ϕ⊤ϕH(\beta, \phi) = -\log p(\beta \mid y) + \tfrac12 \phi^\top \phi

    Give the coefficients a random momentum and follow the gradient of the log posterior with 20 leapfrog steps. Proposals travel far and are still accepted almost every time.

    Watch the trajectories

  4. Expectation propagation

    04

    approximation

    g(β)∝∏igi(β)gi←proj⁡ ⁣[g−i p(yi∣β)]/g−i\begin{aligned} g(\beta) &\propto \textstyle\prod_{i} g_i(\beta) \\ g_i &\leftarrow \operatorname{proj}\!\big[g_{-i}\, p(y_i \mid \beta)\big] / g_{-i} \end{aligned}

    Replace each observation's likelihood by a Gaussian site and refine the sites one at a time by matching the moments of a one-dimensional tilted distribution.

    Step through the sites

What came out

All four roads arrive at the same place

Run on the data as the assignment describes it, the four methods agree on the dose slope to within 0.02. Each extra unit of dose multiplies the odds of improvement by about 1.35, and every method puts essentially all of its posterior mass above zero.

The differences that remain are the instructive ones: the Gaussian approximations are symmetric by construction, while the samplers pick up the slight right skew of the true posterior. HMC needs a fraction of the iterations that random-walk Metropolis does for the same precision.

Posterior for the dose slope β₁, corrected data, values from the original R functions (seed 90125). HMC’s ESS exceeds its 9,000 draws because successive draws are negatively correlated; coda reports it as is.
Methodmean95% interval
Laplace0.302[0.048, 0.555]
MH0.317[0.068, 0.586]
HMC0.315[0.061, 0.580]
EP0.314[0.056, 0.573]

The browser port reproduces every number in this table to at least eleven significant figures; the unit tests compare them with a fresh R run.

Errata

The 2023 submission did not get there

Re-running the original code surfaced data-encoding and likelihood bugs. The algorithms were fine; each one was simply handed a different, wrong version of the data. The explorer keeps that exact behaviour as an “As submitted” mode, and the errata page walks through every bug.

Read the errata Open the 2023 version
Metropolis acceptance
5%
proposals were scaled from the double-counted GLM, about 10× wider than the posterior of the slope being sampled (corrected: 36%)
Smallest effective sample size
1.2
out of 50,000 Metropolis draws: five of seven coefficients never entered the likelihood, so they wandered freely
Improvements HMC and EP saw
7 of 70
the 0/1 outcome label was read as the number improved; the data has 37 of 70

Added in 2026

Can the posterior be trusted?

The 2023 analysis ran one chain per sampler and stopped at a plot. The diagnostics page adds the checks a careful analysis needs, and shows which of them would have caught the 2023 bugs. Every posterior mean and predictive-check p-value carries its Monte Carlo error, and the diagnostics are verified against R’s posterior package.

Open the diagnostics Methods and decision records

  • Four chains per sampler

    Rank-normalised R-hat, bulk and tail ESS, MCSE, trace and rank plots.

  • Posterior predictive checks

    Replicated trials against the observed 14-cell table.

  • Prior sensitivity

    The exact posterior under priors from very tight to flat.

  • Four roads, one exact answer

    Every method against the posterior computed by quadrature.

Under the hood

Faithful by construction

Matching R “roughly” was not the goal. To replay the original chains exactly, the TypeScript port reimplements the pieces of R that the samplers touch, then checks them against a script that runs the submitted R functions.

  • set.seed() and runif()

    R's Mersenne-Twister with its seed scrambling

  • rnorm() and rmvnorm()

    inversion sampling (AS241) and mvtnorm's eigen square root

  • dbinom()

    Loader's saddle-point algorithm, as in R's nmath

  • integrate()

    QUADPACK DQAGS, translated routine by routine

  • glm()

    IRLS with R's starting values and convergence rule

  • effectiveSize(), density()

    coda's AR spectral estimate and R's binned KDE

About this project

Coursework, revived

Subject
MAST90125 Bayesian Statistical Learning
University
The University of Melbourne
When
2023, Semester 2
Assessment
Assignment 3 (individual)
Author
Sunchuangyu (Rin) Huang · sampler and EP code adapted from the MAST90125 lecture notes, as credited in the original
Original stack
R, R Markdown and knitr, mvtnorm, coda
Revived stack
Next.js 16, React 19, TypeScript, Tailwind CSS 4, Web Workers, KaTeX, Vitest; hand-drawn SVG charts

The original submission is preserved unchanged in the repository’s coursework/ folder for reference and academic-integrity transparency. The task is paraphrased here; no assignment specification, lecture material or course dataset files are hosted on this site.

View the source on GitHub