3.2 Sampling to summarize

A posterior distribution contains more information than a single point estimate. In this section, we learn how to use samples to compute credible intervals, posterior means, and other useful summaries.

3.2.1 Intervals of defined mass

A credible interval is an interval that contains a specified fraction of the posterior probability. For example, a 90% credible interval contains 90% of the posterior mass. There are two common types:

Percentile intervals: Use the quantiles of the posterior. A 90% percentile interval goes from the 5th to the 95th percentile.
Highest posterior density (HPD) intervals: The narrowest interval that contains the specified mass. Useful for skewed posteriors.

$\text{Pr}(\theta_{\text{low}} < \theta < \theta_{\text{high}} \mid D) = 0.90$

# Compute credible intervals from posterior samples
# Using the globe tossing posterior (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)

# Percentile intervals

pi50 <- quantile(samples, c(0.25, 0.75))
pi90 <- quantile(samples, c(0.05, 0.95))
pi95 <- quantile(samples, c(0.025, 0.975))

cat("50% percentile interval: [", round(pi50[1], 4), ",", round(pi50[2], 4), "]\n")
cat("90% percentile interval: [", round(pi90[1], 4), ",", round(pi90[2], 4), "]\n")
cat("95% percentile interval: [", round(pi95[1], 4), ",", round(pi95[2], 4), "]\n")

# Plot the intervals

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "p", ylab = "Posterior density",
main = "Credible intervals for p")
abline(v = pi50, col = "firebrick", lty = 2, lwd = 2)
abline(v = pi90, col = "darkgreen", lty = 3, lwd = 2)
abline(v = pi95, col = "purple", lty = 4, lwd = 2)
legend("topright",
legend = c("50% PI", "90% PI", "95% PI"),
col = c("firebrick", "darkgreen", "purple"), lty = c(2, 3, 4), lwd = 2)

Click Run to execute this code in your browser with WebR

3.2.2 Highest posterior density intervals

The highest posterior density (HPD) interval is the narrowest interval that contains the specified probability mass. For symmetric posteriors, it coincides with the percentile interval. For skewed posteriors, it can be different.

# Compare percentile interval vs HPD interval
# For a skewed posterior, they differ

# Build a skewed posterior using a step prior

p*grid <- seq(0, 1, length.out = 1000)
prior <- ifelse(p_grid < 0.5, 0, 1) # step prior: only p >= 0.5
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)

# Percentile interval

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

# HPD interval (simplified): find narrowest interval with 90% mass

# Sort samples and find the narrowest window containing 90%

sorted_samples <- sort(samples)
n <- length(sorted_samples)
window_size <- floor(n * 0.9)
widths <- sorted_samples[(window_size + 1):n] - sorted_samples[1:(n - window_size)]
hpd_idx <- which.min(widths)
hpd90 <- c(sorted_samples[hpd_idx], sorted_samples[hpd_idx + window_size])
cat("90% HPD interval: [", round(hpd90[1], 4), ",", round(hpd90[2], 4), "]\n")
cat("PI width:", round(pi90[2] - pi90[1], 4), "\n")
cat("HPD width:", round(hpd90[2] - hpd90[1], 4), "\n")

# Plot

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "p", ylab = "Posterior density",
main = "Percentile vs HPD interval (skewed posterior)")
abline(v = pi90, col = "firebrick", lty = 2, lwd = 2)
abline(v = hpd90, col = "darkgreen", lty = 3, lwd = 2)
legend("topright",
legend = c("90% PI", "90% HPD"),
col = c("firebrick", "darkgreen"), lty = c(2, 3), lwd = 2)

Click Run to execute this code in your browser with WebR

For this skewed posterior, the HPD interval is narrower than the percentile interval because it concentrates on the region of highest density.

3.2.3 Point estimates

There are three common point estimates from a posterior:

Posterior mean: $E[\theta \mid D] = \frac{1}{S} \sum \theta_s$ . Minimizes expected squared error.
Posterior median: The 50th percentile. Minimizes expected absolute error.
Posterior mode (MAP): The value with highest posterior density. Maximizes posterior probability.

# Compare point estimates for the globe tossing posterior
# Using a symmetric (uniform prior) posterior

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)

post_mean <- mean(samples)
post_median <- median(samples)
post_mode <- p_grid[which.max(posterior)]

cat("Posterior mean:", round(post_mean, 4), "\n")
cat("Posterior median:", round(post_median, 4), "\n")
cat("Posterior mode (MAP):", round(post_mode, 4), "\n")
cat("MLE (w/n):", round(6/9, 4), "\n")

# Plot

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "p", ylab = "Posterior density",
main = "Point estimates from the posterior")
abline(v = post_mean, col = "firebrick", lty = 2, lwd = 2)
abline(v = post_median, col = "darkgreen", lty = 3, lwd = 2)
abline(v = post_mode, col = "purple", lty = 4, lwd = 2)
legend("topright",
legend = c(paste("Mean:", round(post_mean, 3)),
paste("Median:", round(post_median, 3)),
paste("Mode:", round(post_mode, 3))),
col = c("firebrick", "darkgreen", "purple"), lty = c(2, 3, 4), lwd = 2)

Click Run to execute this code in your browser with WebR

For a symmetric posterior, all three estimates are similar. For a skewed posterior, they differ: the mean is pulled toward the tail, the median is in the middle, and the mode is at the peak.

3.2.4 Posterior probabilities

We can compute the probability of any event by counting the fraction of samples that satisfy it:

$P(\theta \in A \mid D) = \frac{1}{S} \sum_{s=1}^S \mathbf{1}[\theta_s \in A]$

# Compute posterior probabilities of various events
# Using the globe tossing posterior

cat("P(p < 0.5) =", round(mean(samples < 0.5), 4), "\n")
cat("P(p > 0.5) =", round(mean(samples > 0.5), 4), "\n")
cat("P(0.5 < p < 0.8) =", round(mean(samples > 0.5 & samples < 0.8), 4), "\n")
cat("P(p > 0.6) =", round(mean(samples > 0.6), 4), "\n")
cat("P(p > 0.7) =", round(mean(samples > 0.7), 4), "\n")

# These are exact posterior probabilities, not p-values

# They tell us exactly what we want to know about p

Click Run to execute this code in your browser with WebR

Rethinking

Intervals versus point estimates. A point estimate is a summary of the posterior, but it discards uncertainty. An interval preserves some uncertainty, but it still reduces the posterior to two numbers. The full posterior distribution is always the most complete answer. When reporting results, always show the posterior distribution (or a histogram of samples) alongside any summaries. The summaries are for communication; the distribution is for understanding.

Summary of Section 3.2

Credible intervals contain a specified fraction of the posterior mass.
Percentile intervals use posterior quantiles; HPD intervals are the narrowest intervals.
Point estimates include the mean (minimizes squared error), median (minimizes absolute error), and mode (MAP).
Posterior probabilities are computed by counting the fraction of samples that satisfy an event.
The full posterior distribution is always the most complete answer; summaries are for communication.

Sampling from a grid-approximate posterior

Sampling to simulate prediction