3.5 Práctica
Los problemas de práctica del Capítulo 3 le piden muestrear de posteriors, calcular intervalos de credibilidad y simular distribuciones predictivas a posteriori. Enfatizan que la posterior es una distribución, no un punto, y que los resúmenes deben preservar la incertidumbre.
Problemas fáciles
3E1. ¿Cuánta probabilidad a posteriori hay por debajo de ?
# 3E1: P(p < 0.2) from the globe tossing 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)
cat("P(p < 0.2) =", round(mean(samples < 0.2), 4), "\n")3E2. ¿Cuánta probabilidad a posteriori hay por encima de ?
3E3. ¿Cuánta probabilidad a posteriori hay entre y ?
3E4. ¿Cuál es el percentil 20 de las muestras posteriores?
# 3E4: 20th percentile of posterior
cat("20th percentile:", round(quantile(samples, 0.2), 4), "\n")3E5. ¿Cuál es el percentil 80 de las muestras posteriores?
Problemas medios
3M1. Suponga que los datos del lanzamiento del globo hubieran sido 8 aguas en 15 lanzamientos. Construya la posterior y calcule el intervalo HPD del 90%.
# 3M1: 8 waters in 15 tosses
p_grid <- seq(0, 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood <- dbinom(8, size = 15, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood * prior)
set.seed(42)
samples <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior)
# 90% percentile interval
pi90 <- quantile(samples, c(0.05, 0.95))
cat("90% percentile interval: [", round(pi90[1], 4), ",", round(pi90[2], 4), "]\n")
# 90% HPD interval
sorted <- sort(samples)
n <- length(sorted)
window <- floor(n * 0.9)
widths <- sorted[(window+1):n] - sorted[1:(n-window)]
hpd_idx <- which.min(widths)
hpd90 <- c(sorted[hpd_idx], sorted[hpd_idx + window])
cat("90% HPD interval: [", round(hpd90[1], 4), ",", round(hpd90[2], 4), "]\n")
# Plot
plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "p", ylab = "Posterior density",
main = "3M1: 8 waters in 15 tosses")
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)3M2. Extraiga 10.000 muestras de la posterior para 8 aguas en 15 lanzamientos. Simule muestras predictivas a posteriori para 15 lanzamientos. ¿Cuál es la probabilidad predicha de observar 8 aguas?
3M3. Compare la distribución predictiva a posteriori con los datos observados. ¿Es 8 aguas un valor típico?
Problemas difíciles
3H1. Usando la posterior de 3M1, calcule la probabilidad a posteriori de que .
3H2. Simule 10.000 conjuntos de datos predictivos a posteriori de 15 lanzamientos cada uno. Calcule la proporción de simulaciones que producen 8 o más aguas. Compare esto con los datos observados.
# 3H2: Posterior predictive check for 8/15
p_grid <- seq(0, 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood <- dbinom(8, size = 15, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood * prior)
set.seed(42)
samples <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior)
# Simulate posterior predictive
pred <- rbinom(1e4, size = 15, prob = samples)
cat("P(simulated >= 8):", round(mean(pred >= 8), 4), "\n")
cat("P(simulated <= 8):", round(mean(pred <= 8), 4), "\n")
cat("Observed: 8 waters in 15 tosses\n")
hist(pred, breaks = -0.5:15.5, col = "steelblue", border = "white",
xlab = "Number of waters", main = "3H2: Posterior predictive (8/15)")
abline(v = 8, col = "firebrick", lwd = 3)3H3. Suponga que observa 3 aguas en 3 lanzamientos. Construya la posterior y simule muestras predictivas a posteriori. ¿Cuál es la probabilidad predicha de observar 3 aguas en 3 lanzamientos futuros? ¿Cómo se compara con el caso con más datos?