3.3 Muestreo para simular predicciones

Las muestras no solo son útiles para resumir la posterior — también son esenciales para hacer predicciones. Una distribución predictiva a posteriori genera nuevos datos consistentes con la posterior, permitiéndonos verificar el modelo y pronosticar observaciones futuras.

3.3.1 La distribución predictiva a posteriori

La distribución predictiva a posteriori es la distribución de nuevos datos $y_\text{new}$ , dados los datos observados $D$ :

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

Esta integral promedia la verosimilitud sobre la posterior, produciendo una distribución que cuenta tanto con la incertidumbre del parámetro como con la variabilidad muestral. Usando muestras, la aproximamos por:

Extraer $\theta_s$ de la posterior.
Extraer $y_\text{new,s}$ de la verosimilitud $p(y_\text{new} \mid \theta_s)$ .
Repetir.

La colección $\{y_\text{new,1}, \ldots, y_\text{new,S}\}$ es una muestra de la distribución predictiva a posteriori.

# 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

La distribución predictiva a posteriori muestra qué nuevos datos esperaríamos ver. Los datos observados (6 aguas) deberían ser un valor típico bajo esta distribución. Si no lo es, el modelo puede ser inadecuado.

3.3.2 Verificación del modelo con predicciones a posteriori

Las verificaciones predictivas a posteriori comparan los datos observados con la distribución predictiva a posteriori. Si el modelo es adecuado, los datos observados deberían verse como una extracción típica de la distribución predictiva.

# 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 Inflación de ceros y desajuste del modelo

Supongamos que observamos 0 aguas en 9 lanzamientos — un valor extremo. ¿Seguiría ajustando nuestro modelo?

# 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 Prediciendo observaciones futuras

También podemos usar la distribución predictiva a posteriori para pronosticar datos futuros. Por ejemplo, si lanzamos el globo 10 veces más, ¿cuántas aguas deberíamos esperar?

# 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

La predicción es la prueba definitiva. Un modelo que ajusta el pasado pero no puede predecir el futuro es inútil. Las verificaciones predictivas a posteriori son la forma más simple de validación de predicciones: verificamos si el modelo puede reproducir los patrones que ya hemos visto. Pero la verdadera prueba es la predicción fuera de muestra — ¿puede el modelo predecir datos que no ha visto? Este es el enfoque del Capítulo 6, donde usamos criterios de información para estimar la precisión predictiva fuera de muestra.

Overthinking

La predictiva a posteriori como distribución de mezcla. La distribución predictiva a posteriori es una distribución de mezcla: es un promedio ponderado de verosimilitudes, donde los pesos son las probabilidades a posteriori:

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

Esto es diferente de sustituir una sola estimación puntual. Si usáramos solo la media a posteriori $\hat{\theta}$ , obtendríamos $p(y_\text{new} \mid \hat{\theta})$ , que ignora la incertidumbre del parámetro. La distribución predictiva a posteriori es más ancha porque integra sobre la posterior, contando tanto con la variabilidad en $\theta$ como con la variabilidad en $y_\text{new}$ dado $\theta$ .

Para el modelo binomial, la predictiva a posteriori es una distribución Beta-Binomial. Si la posterior es $\text{Beta}(\alpha, \beta)$ , entonces la distribución predictiva para $y_\text{new}$ éxitos en $n$ ensayos es:

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

Esta distribución tiene colas más pesadas que la binomial, reflejando la incertidumbre adicional de no conocer $p$ exactamente.

Resumen de la Sección 3.3

La distribución predictiva a posteriori $p(y_\text{new} \mid D)$ promedia la verosimilitud sobre la posterior.
Para muestrear de ella: extraiga $\theta_s$ de la posterior, luego extraiga $y_\text{new}$ de la verosimilitud.
Las verificaciones predictivas a posteriori comparan datos observados con datos simulados para evaluar el ajuste del modelo.
La predictiva a posteriori es una distribución de mezcla — más ancha que sustituir una estimación puntual.
La predicción es la prueba definitiva de la utilidad de un modelo.

Muestreo para resumir

Resumen