3.2 Muestreo para resumir

Una distribución a posteriori contiene más información que una sola estimación puntual. En esta sección, aprendemos a usar muestras para calcular intervalos de credibilidad, medias a posteriori y otros resúmenes útiles.

3.2.1 Intervalos de masa definida

Un intervalo de credibilidad es un intervalo que contiene una fracción especificada de la probabilidad a posteriori. Por ejemplo, un intervalo de credibilidad del 90% contiene el 90% de la masa a posteriori. Hay dos tipos comunes:

Intervalos de percentiles: Usan los cuantiles de la posterior. Un intervalo de percentil del 90% va del percentil 5 al 95.
Intervalos de densidad a posteriori más alta (HPD): El intervalo más estrecho que contiene la masa especificada. Útiles para posteriors asimétricas.

$\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 Intervalos de densidad a posteriori más alta

El intervalo de densidad a posteriori más alta (HPD) es el intervalo más estrecho que contiene la masa de probabilidad especificada. Para posteriors simétricas, coincide con el intervalo de percentiles. Para posteriors asimétricas, puede ser diferente.

# 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

Para esta posterior asimétrica, el intervalo HPD es más estrecho que el intervalo de percentiles porque se concentra en la región de mayor densidad.

3.2.3 Estimaciones puntuales

Hay tres estimaciones puntuales comunes de una posterior:

Media a posteriori: $E[\theta \mid D] = \frac{1}{S} \sum \theta_s$ . Minimiza el error cuadrático esperado.
Mediana a posteriori: El percentil 50. Minimiza el error absoluto esperado.
Moda a posteriori (MAP): El valor con mayor densidad a posteriori. Maximiza la probabilidad a posteriori.

# 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

Para una posterior simétrica, las tres estimaciones son similares. Para una posterior asimétrica, difieren: la media es atraída hacia la cola, la mediana está en el medio, y la moda está en el pico.

3.2.4 Probabilidades a posteriori

Podemos calcular la probabilidad de cualquier evento contando la fracción de muestras que lo satisfacen:

$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

Intervalos versus estimaciones puntuales. Una estimación puntual es un resumen de la posterior, pero descarta la incertidumbre. Un intervalo preserva algo de incertidumbre, pero aún reduce la posterior a dos números. La distribución a posteriori completa es siempre la respuesta más completa. Al reportar resultados, siempre muestre la distribución a posteriori (o un histograma de muestras) junto con cualquier resumen. Los resúmenes son para comunicación; la distribución es para entendimiento.

Resumen de la Sección 3.2

Los intervalos de credibilidad contienen una fracción especificada de la masa a posteriori.
Los intervalos de percentiles usan cuantiles posteriores; los intervalos HPD son los intervalos más estrechos.
Las estimaciones puntuales incluyen la media (minimiza error cuadrático), la mediana (minimiza error absoluto) y la moda (MAP).
Las probabilidades a posteriori se calculan contando la fracción de muestras que satisfacen un evento.
La distribución a posteriori completa es siempre la respuesta más completa; los resúmenes son para comunicación.

Muestreo desde una posterior aproximada por cuadrícula

Muestreo para simular predicciones