2.4 Making the model go

Once a model is specified — likelihood, parameters, and prior — we need to compute the posterior. This is the practical heart of Bayesian statistics. In this section, we explore three strategies for computing the posterior: grid approximation, quadratic approximation, and Markov Chain Monte Carlo.

2.4.1 Grid approximation

Grid approximation is the simplest method. We lay out a grid of possible parameter values, compute the posterior at each grid point, and normalize. This is exactly what we did in Section 2.1.

The advantage of grid approximation is that it is exact (up to the grid resolution) and transparent — we can see every step. The disadvantage is that it scales poorly: with $k$ parameters and $g$ grid points per parameter, we need $g^k$ computations. For a model with 10 parameters and 100 grid points each, that is $100^{10} = 10^{20}$ computations — infeasible.

# Grid approximation for the globe tossing problem
# Data: 6 waters in 9 tosses

# 1. Define the grid

p_grid <- seq(from = 0, to = 1, length.out = 100)

# 2. Define the prior (uniform)

prior <- rep(1, 100)

# 3. Compute the likelihood

likelihood <- dbinom(6, size = 9, prob = p_grid)

# 4. Compute the unstandardized posterior

unstd_posterior <- likelihood * prior

# 5. Standardize

posterior <- unstd_posterior / sum(unstd_posterior)

# 6. Plot

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "Probability of water (p)",
ylab = "Posterior probability",
main = "Grid approximation: 100 points")
polygon(c(p_grid, rev(p_grid)),
c(posterior, rep(0, length(posterior))),
col = rgb(0.27, 0.51, 0.71, 0.3), border = NA)

# 7. Summarize

cat("Posterior mean:", round(sum(p_grid * posterior), 4), "\n")
cat("Posterior mode:", round(p_grid[which.max(posterior)], 4), "\n")

Click Run to execute this code in your browser with WebR

Grid approximation works well for models with 1-2 parameters. Beyond that, we need more efficient methods.

2.4.2 Quadratic approximation

Quadratic approximation (also called Laplace approximation) approximates the posterior with a multivariate Gaussian distribution. The idea is to find the peak of the posterior (the maximum a posteriori or MAP estimate) and then approximate the posterior as a Gaussian centered at the peak, with covariance determined by the curvature of the log-posterior.

The log-posterior near the peak is approximately quadratic:

$\log p(\theta \mid D) \approx \log p(\hat{\theta} \mid D) - \frac{1}{2}(\theta - \hat{\theta})^T H (\theta - \hat{\theta})$

where $\hat{\theta}$ is the MAP estimate and $H$ is the Hessian matrix (matrix of second derivatives) of the log-posterior. The approximate posterior is:

$p(\theta \mid D) \approx \text{Normal}(\hat{\theta}, H^{-1})$

# Quadratic approximation for the globe tossing problem
# Compare with grid approximation

# Grid approximation (reference)

p*grid <- seq(0, 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood <- dbinom(6, size = 9, prob = p_grid)
post_grid <- (likelihood * prior) / sum(likelihood _ prior)

# Quadratic approximation

# The MAP estimate for a binomial with uniform prior is w/n

w <- 6; n <- 9
p_map <- w / n

# The standard deviation of the quadratic approximation

# For binomial: sigma = sqrt(p*(1-p)/n)

sigma_quad <- sqrt(p_map * (1 - p_map) / n)

cat("MAP estimate:", round(p_map, 4), "\n")
cat("Quadratic approximation sigma:", round(sigma_quad, 4), "\n")

# Compare grid and quadratic

plot(p_grid, post_grid * 1000, type = "l", lwd = 2, col = "steelblue",
xlab = "Probability of water (p)",
ylab = "Density",
main = "Grid vs Quadratic approximation")
curve(dnorm(x, mean = p_map, sd = sigma_quad), from = 0, to = 1,
add = TRUE, lwd = 2, col = "firebrick", lty = 2)
legend("topright",
legend = c("Grid approximation", "Quadratic approximation"),
col = c("steelblue", "firebrick"), lty = c(1, 2), lwd = 2)

Click Run to execute this code in your browser with WebR

The quadratic approximation is very close to the grid approximation when the posterior is roughly bell-shaped. The advantage of the quadratic approximation is that it is fast — we only need to find the peak and compute the curvature. The disadvantage is that it fails when the posterior is skewed, multimodal, or has bounded support (like $p \in [0, 1]$ ).

2.4.3 Markov Chain Monte Carlo

Markov Chain Monte Carlo (MCMC) is the most flexible method. Instead of computing the posterior on a grid or approximating it with a Gaussian, MCMC draws samples from the posterior. These samples can be used to compute any quantity of interest (means, intervals, predictions) by taking averages.

The key idea is to construct a Markov chain — a sequence of parameter values where each value depends only on the previous one — whose stationary distribution is the posterior. After the chain has converged, the samples are distributed according to the posterior.

MCMC scales to models with hundreds or thousands of parameters, making it the workhorse of modern Bayesian statistics. We will study MCMC in detail in Chapter 8.

# MCMC from scratch: Metropolis algorithm for the globe tossing problem
# Data: 6 waters in 9 tosses

set.seed(42)
n_samples <- 1000
p_current <- 0.5 # starting value
samples <- numeric(n_samples)

for (i in 1:n_samples) {

# Propose a new value by adding a small random jump

p_proposal <- p_current + runif(1, -0.1, 0.1)

# Ensure the proposal stays in [0, 1]

if (p_proposal < 0) p_proposal <- abs(p_proposal)
if (p_proposal > 1) p_proposal <- 2 - p_proposal

# Compute the acceptance ratio

# Posterior is proportional to likelihood * prior

# Prior is uniform, so it cancels

log_post_current <- dbinom(6, 9, p_current, log = TRUE)
log_post_proposal <- dbinom(6, 9, p_proposal, log = TRUE)

# Accept or reject

if (log(runif(1)) < (log_post_proposal - log_post_current)) {
p_current <- p_proposal
}
samples[i] <- p_current
}

# Plot the samples

par(mfrow = c(1, 2))

# Trace plot

plot(1:n_samples, samples, type = "l", col = "steelblue",
xlab = "Iteration", ylab = "p",
main = "MCMC trace plot")

# Histogram of samples

hist(samples, breaks = 30, col = "steelblue", border = "white",
freq = FALSE, xlab = "p",
main = "MCMC posterior samples")

# Overlay grid approximation

p_grid <- seq(0, 1, length.out = 100)
post_grid <- dbinom(6, 9, p_grid) / sum(dbinom(6, 9, p_grid))
lines(p_grid, post_grid * 100, lwd = 2, col = "firebrick")

par(mfrow = c(1, 1))

cat("MCMC posterior mean:", round(mean(samples), 4), "\n")
cat("MCMC posterior SD:", round(sd(samples), 4), "\n")
cat("MCMC 95% interval: [", round(quantile(samples, 0.025), 4), ",",
round(quantile(samples, 0.975), 4), "]\n")

Click Run to execute this code in your browser with WebR

The trace plot shows the chain exploring the parameter space. The histogram of samples approximates the posterior distribution. With more samples, the approximation improves. MCMC is the most powerful method, but it requires careful diagnosis (Chapter 8).

2.4.4 Comparing the three methods

| Method | Exact? | Scalable? | Fast? | When to use | | --------- | ---------------- | ---------- | ---------------- | ------------------------- | | Grid | Yes (up to grid) | No ( $g^k$ ) | Slow for $k > 2$ | Teaching, 1-2 parameters | | Quadratic | Approximate | Yes | Very fast | Posteriors near Gaussian | | MCMC | Asymptotically | Yes | Moderate | Complex models, any shape |

Rethinking

The small world and the large world. All three methods operate in the small world — the mathematical model we have built. The small world is self-consistent: given the model, the posterior is well-defined, and the methods approximate it. But the large world — reality — may not match the model. A perfect computation of the wrong posterior is still wrong. The methods in this section tell us how to compute the posterior accurately. They do not tell us whether the model is correct. That is the job of model criticism (Chapter 6) and scientific judgment.

Overthinking

The Hessian and the quadratic approximation. The Hessian matrix $H$ of the log-posterior contains the second derivatives:

$H_{ij} = \frac{\partial^2 \log p(\theta \mid D)}{\partial \theta_i \partial \theta_j}$

The inverse of the Hessian $H^{-1}$ is the covariance matrix of the quadratic approximation. For a model with a single parameter, $H$ is a scalar and $\sigma = 1/\sqrt{-H}$ . For the binomial model with $w$ successes in $n$ trials and a uniform prior:

$\log p(p \mid D) = w \log p + (n-w) \log(1-p) + \text{const}$

$\frac{d^2 \log p}{dp^2} = -\frac{w}{p^2} - \frac{n-w}{(1-p)^2}$

At the MAP $\hat{p} = w/n$ :

$H = -\frac{w}{(w/n)^2} - \frac{n-w}{(1-w/n)^2} = -\frac{n^2}{w} - \frac{n^2}{n-w} = -\frac{n^3}{w(n-w)}$

$\sigma = \frac{1}{\sqrt{n^3 / (w(n-w))}} = \sqrt{\frac{w(n-w)}{n^3}} = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$

This is the familiar standard error of a proportion, confirming that the quadratic approximation recovers the classical result.

Summary of Section 2.4

Grid approximation computes the posterior at discrete parameter values. Exact but scales as $g^k$ .
Quadratic approximation fits a Gaussian to the posterior peak. Fast and accurate for bell-shaped posteriors.
MCMC draws samples from the posterior. Flexible and scales to complex models, but requires diagnosis.
The three methods trade off accuracy, speed, and scalability.
All methods operate in the small world of the model — they do not guarantee the model is correct.

Components of the model

Summary