MCMC and the Metropolis Sampler

Lesson 2 ended on a confession: conjugate updating is beautiful, and most real models do not qualify. One predictor, one unknown spread, one hierarchy, and there is no family to add your way into and no formula for the posterior at all. This lesson builds the escape route: an algorithm from 1953 that draws samples from a posterior it can never write down, and that still powers every serious Bayesian tool today. You will build it yourself, in fifteen lines of base R.

By the end of this lesson you will be able to:

• Say precisely which part of Bayes' theorem becomes impossible in real models, and why drawing samples gets around it
• Turn a bag of posterior draws into every deliverable Lesson 2 produced: a mean, a credible interval, any probability
• Build the Metropolis sampler from scratch and validate it against an answer you know exactly
• Read a trace plot, discard burn-in, and tune the proposal width by watching the acceptance rate
• Point the same fifteen lines at a model with no closed form, by changing one function

Prerequisites: Lessons 1 and 2 of this course (prior, likelihood, posterior, Asha's Beta(16, 104) conversion posterior, the Normal-Normal update for her average order value), plus base R functions and for loops.

Below is where we are headed: a random walk exploring a posterior it cannot see, leaving behind a histogram of everywhere it stood. Toggle the step size and watch the walk crawl, roam or stall. By the end of this lesson you will have built this machine from nothing.

Where the formulas run out

The day after the redesign launch, Asha's checkout page takes 10 orders. Nine are ordinary baskets, between $39 and $71. One is a $237 corporate bulk buy, a gift shipment for someone's whole office. Her order-value model from Lesson 2 assumed every order is a draw from one normal curve with an $18 spread. On that assumption a $237 order is a ten-standard-deviations miracle, so the model will wrench its estimate of the average toward it. What Asha wants instead is a likelihood that treats big orders as rare but unremarkable: a bell curve with heavier tails. Statisticians keep exactly that curve on the shelf. It is called the Student t.

And right there, the machinery of Lessons 1 and 2 dies. Write the update anyway, for \(\mu\), the true long-run average order value:

\[ p(\mu \mid \text{data}) \;=\; \frac{p(\text{data} \mid \mu)\; p(\mu)}{\displaystyle\int p(\text{data} \mid \mu')\, p(\mu')\, d\mu'} \]

The numerator is the old friend: likelihood times prior, a number you can compute at any single candidate \(\mu\) in a microsecond. The denominator is that same numerator ADDED UP across every value the average could possibly be (\(\mu'\) just means "each candidate in turn", and the integral sign means "sum over all of them"). For a matched prior-likelihood pair, algebra does that infinite sum for you: that was conjugacy. A normal prior with a t likelihood is not a matched pair. Nobody has a formula for its total. Not "the formula is hard": it does not exist in closed form.

Could the grid from Lesson 1 save us? For one unknown, honestly, yes. But real models stack unknowns, and a grid grows exponentially: