3.3 Sampling to simulate prediction

Samples are not only useful for summarizing the posterior — they are also essential for making predictions. A posterior predictive distribution generates new data consistent with the posterior, allowing us to check the model and forecast future observations.

3.3.1 The posterior predictive distribution

The posterior predictive distribution is the distribution of new data $y_\text{new}$ , given the observed data $D$ :

$p(y_\text{new} \mid D) = \int p(y_\text{new} \mid \theta) \cdot p(\theta \mid D) \, d\theta$

This integral averages the likelihood over the posterior, producing a distribution that accounts for both parameter uncertainty and sampling variability. Using samples, we approximate it by:

Draw $\theta_s$ from the posterior.
Draw $y_\text{new,s}$ from the likelihood $p(y_\text{new} \mid \theta_s)$ .
Repeat.

The collection $\{y_\text{new,1}, \ldots, y_\text{new,S}\}$ is a sample from the posterior predictive distribution.

# Posterior predictive simulation for the globe tossing problem
# Data: 6 waters in 9 tosses

# Build the posterior and sample

p*grid <- seq(0, 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood <- dbinom(6, size = 9, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood _ prior)
set.seed(42)
samples <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior)

# Posterior predictive: simulate new data

# For each posterior sample, simulate 9 tosses

posterior_predictive <- rbinom(1e4, size = 9, prob = samples)

# Tabulate the results

pred_table <- table(posterior_predictive) / length(posterior_predictive)
cat("Posterior predictive distribution:\n")
print(round(pred_table, 4))

# Compare with observed data (6 waters)

cat("\nObserved: 6 waters\n")
cat("P(6 waters | model):", round(pred_table["6"], 4), "\n")

# Plot

barplot(pred_table, col = "steelblue", border = "white",
xlab = "Number of waters in 9 tosses",
ylab = "Predictive probability",
main = "Posterior predictive distribution")
abline(v = which(names(pred_table) == "6") + 0.7,
lty = 2, col = "firebrick", lwd = 2)
legend("topright", legend = "Observed (6)", col = "firebrick", lty = 2, lwd = 2)

Click Run to execute this code in your browser with WebR

The posterior predictive distribution shows what new data we would expect to see. The observed data (6 waters) should be a typical value under this distribution. If it is not, the model may be inadequate.

3.3.2 Model checking with posterior predictions

Posterior predictive checks compare the observed data to the posterior predictive distribution. If the model is adequate, the observed data should look like a typical draw from the predictive distribution.

# Posterior predictive check: is 6 waters typical?
# Compare observed data to simulated data

# Simulate many datasets from the posterior predictive

n_sim <- 1000
simulated_datasets <- rbinom(n_sim, size = 9, prob = samples[1:n_sim])

# Summary statistics of simulated data

cat("Simulated data summary:\n")
cat("Mean:", round(mean(simulated_datasets), 3), "\n")
cat("SD:", round(sd(simulated_datasets), 3), "\n")
cat("Observed: 6\n")
cat("P(simulated >= observed):", round(mean(simulated_datasets >= 6), 4), "\n")
cat("P(simulated <= observed):", round(mean(simulated_datasets <= 6), 4), "\n")

# If both probabilities are reasonable (not close to 0 or 1),

# the model fits the data

# Visualize

hist(simulated_datasets, breaks = -0.5:9.5, col = "steelblue",
border = "white", xlab = "Number of waters",
main = "Posterior predictive check")
abline(v = 6, col = "firebrick", lwd = 3)
legend("topright", legend = "Observed (6)", col = "firebrick", lwd = 3)

Click Run to execute this code in your browser with WebR

3.3.3 Zero-inflation and model misfit

Suppose we observed 0 waters in 9 tosses — an extreme value. Would our model still fit?

# What if we observed 0 waters instead of 6?
# Would the model still fit?

# Posterior for 0 waters in 9 tosses

p*grid <- seq(0, 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood_0 <- dbinom(0, size = 9, prob = p_grid)
posterior_0 <- (likelihood_0 * prior) / sum(likelihood*0 * prior)
set.seed(42)
samples_0 <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior_0)

# Posterior predictive

pred_0 <- rbinom(1e4, size = 9, prob = samples_0)

cat("If we observed 0 waters:\n")
cat("Posterior mean of p:", round(mean(samples_0), 4), "\n")
cat("Posterior predictive mean:", round(mean(pred_0), 4), "\n")
cat("P(0 waters in new data):", round(mean(pred_0 == 0), 4), "\n")

# The model can accommodate 0 waters by pushing p toward 0

# But this might indicate a problem: maybe the globe never lands on water

# This is where model criticism becomes important

Click Run to execute this code in your browser with WebR

3.3.4 Predicting future observations

We can also use the posterior predictive distribution to forecast future data. For example, if we toss the globe 10 more times, how many waters should we expect?

# Predict future observations: 10 more tosses
# Using the posterior from 6 waters in 9 tosses

p*grid <- seq(0, 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood <- dbinom(6, size = 9, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood _ prior)
set.seed(42)
samples <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior)

# Predict 10 future tosses

future_tosses <- rbinom(1e4, size = 10, prob = samples)

cat("Prediction for 10 future tosses:\n")
cat("Expected number of waters:", round(mean(future_tosses), 2), "\n")
cat("90% prediction interval: [",
round(quantile(future_tosses, 0.05), 0), ",",
round(quantile(future_tosses, 0.95), 0), "]\n")

# Plot

hist(future_tosses, breaks = -0.5:10.5, col = "steelblue",
border = "white", xlab = "Number of waters in 10 future tosses",
main = "Posterior predictive: 10 future tosses")

Click Run to execute this code in your browser with WebR

Rethinking

Prediction is the ultimate test. A model that fits the past but cannot predict the future is useless. Posterior predictive checks are the simplest form of prediction validation: we check whether the model can reproduce the patterns we have already seen. But the real test is out-of-sample prediction — can the model predict data it has not seen? This is the focus of Chapter 6, where we use information criteria to estimate out-of-sample predictive accuracy.

Overthinking

The posterior predictive as a mixture distribution. The posterior predictive distribution is a mixture distribution: it is a weighted average of likelihoods, where the weights are the posterior probabilities:

$p(y_\text{new} \mid D) = E_{\theta \mid D}[p(y_\text{new} \mid \theta)]$

This is different from plugging in a single point estimate. If we used only the posterior mean $\hat{\theta}$ , we would get $p(y_\text{new} \mid \hat{\theta})$ , which ignores parameter uncertainty. The posterior predictive distribution is wider because it integrates over the posterior, accounting for both the variability in $\theta$ and the variability in $y_\text{new}$ given $\theta$ .

For the binomial model, the posterior predictive is a Beta-Binomial distribution. If the posterior is $\text{Beta}(\alpha, \beta)$ , then the predictive distribution for $y_\text{new}$ successes in $n$ trials is:

$y_\text{new} \sim \text{Beta-Binomial}(n, \alpha, \beta)$

This distribution has heavier tails than the binomial, reflecting the extra uncertainty from not knowing $p$ exactly.

Summary of Section 3.3

The posterior predictive distribution $p(y_\text{new} \mid D)$ averages the likelihood over the posterior.
To sample from it: draw $\theta_s$ from the posterior, then draw $y_\text{new}$ from the likelihood.
Posterior predictive checks compare observed data to simulated data to assess model fit.
The posterior predictive is a mixture distribution — wider than plugging in a point estimate.
Prediction is the ultimate test of a model's usefulness.

Sampling to summarize

Summary