3.1 Sampling from a grid-approximate posterior

Sampling from the posterior is one of the most useful techniques in Bayesian statistics. Once we have a sample of parameter values, we can approximate any posterior quantity by computing the corresponding function on the sample. This turns hard integrals into simple averages.

3.1.1 Why sample?

The posterior distribution is a probability distribution over the parameters. In principle, we could compute any quantity we want by integration:

$E[f(\theta)] = \int f(\theta) \, p(\theta \mid D) \, d\theta$

But in practice, this integral is often impossible to compute analytically. Sampling solves this problem. If we have $S$ samples $\theta_1, \theta_2, \ldots, \theta_S$ from the posterior, we can approximate the integral by the sample average:

$\hat{E}[f(\theta)] \approx \frac{1}{S} \sum_{s=1}^{S} f(\theta_s)$

This is the Monte Carlo principle: replace integration with averaging. The law of large numbers guarantees that the approximation improves as $S$ increases.

3.1.2 Sampling from a grid-approximate posterior

In Chapter 2, we computed the posterior on a grid of parameter values. Each grid point has a posterior probability. To sample from this distribution, we simply draw grid points with probabilities proportional to their posterior weights.

# Sample from a grid-approximate posterior
# Globe tossing: 6 waters in 9 tosses

# 1. Build the grid and posterior

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

# 2. Sample 10,000 values from the grid

set.seed(42)
samples <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior)

# 3. Plot the samples

par(mfrow = c(1, 2))

# Grid posterior

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "p", ylab = "Posterior density",
main = "Grid posterior")

# Sampled posterior

hist(samples, breaks = 50, col = "steelblue", border = "white",
freq = FALSE, xlab = "p",
main = "Sampled posterior (10,000 draws)")
lines(p_grid, posterior * 1000, lwd = 2, col = "firebrick")

par(mfrow = c(1, 1))

cat("Number of samples:", length(samples), "\n")
cat("Sample mean:", round(mean(samples), 4), "\n")
cat("Sample SD:", round(sd(samples), 4), "\n")

Click Run to execute this code in your browser with WebR

The sample function draws 10,000 values from the grid, with probabilities proportional to the posterior. Each draw is a plausible value of $p$ . The histogram of these draws approximates the posterior distribution.

3.1.3 What can we do with samples?

Once we have samples from the posterior, we can:

Estimate the posterior mean: mean(samples)
Estimate credible intervals: quantile(samples, c(0.05, 0.95))
Estimate posterior probabilities: mean(samples > 0.5)
Simulate predictions: rbinom(1e4, size = 9, prob = samples)
Compute any function of parameters: mean(samples^2) for $E[p^2]$

The key insight is that every question about the posterior becomes an empirical question about the samples. We never need to compute another integral — we just compute on the samples and take the average.

# Everything we can do with posterior samples
# Using the samples from the previous code block

# Posterior mean

cat("Posterior mean of p:", round(mean(samples), 4), "\n")

# Posterior median

cat("Posterior median:", round(median(samples), 4), "\n")

# 90% credible interval (percentile)

ci90 <- quantile(samples, c(0.05, 0.95))
cat("90% percentile interval: [", round(ci90[1], 4), ",", round(ci90[2], 4), "]\n")

# Posterior probability that p > 0.5

cat("P(p > 0.5):", round(mean(samples > 0.5), 4), "\n")

# Posterior probability that p > 0.6

cat("P(p > 0.6):", round(mean(samples > 0.6), 4), "\n")

# Expected value of p^2

cat("E[p^2]:", round(mean(samples^2), 4), "\n")

# Variance of p

cat("Var(p):", round(var(samples), 6), "\n")

Click Run to execute this code in your browser with WebR

Rethinking

Sampling is the universal solvent. The beauty of sampling is that it unifies all posterior computations. Whether you want a mean, an interval, a probability, or a prediction, the answer is the same: compute the function on the samples and take the average. This is why MCMC (Chapter 8) is so powerful — it produces samples from any posterior, no matter how complex. Once you have samples, the rest is arithmetic.

Overthinking

Why does sampling work? The law of large numbers guarantees that if $\theta_1, \ldots, \theta_S$ are i.i.d. samples from $p(\theta \mid D)$ , then:

$\frac{1}{S} \sum_{s=1}^S f(\theta_s) \xrightarrow{a.s.} E[f(\theta)] = \int f(\theta) \, p(\theta \mid D) \, d\theta$

The rate of convergence is $O(1/\sqrt{S})$ , by the central limit theorem. This means that to halve the Monte Carlo error, we need four times as many samples. With 10,000 samples, the error is typically about 1% of the standard deviation — more than enough for practical purposes.

Importantly, the convergence rate does not depend on the dimension of $\theta$ . This is why sampling scales to high-dimensional models while grid approximation does not. Grid approximation requires $g^k$ points for $k$ parameters, but sampling requires the same number of samples regardless of $k$ (though the mixing of the Markov chain may slow down in high dimensions).

Summary of Section 3.1

Sampling replaces integration with averaging: $\hat{E}[f(\theta)] \approx \frac{1}{S} \sum f(\theta_s)$ .
To sample from a grid posterior, draw grid points with probability proportional to their posterior weight.
Once we have samples, every posterior question becomes an empirical question.
The law of large numbers guarantees convergence; the rate is $O(1/\sqrt{S})$ .
Sampling scales to high dimensions, unlike grid approximation.

Practice

Sampling to summarize