3.1 Muestreo desde una posterior aproximada por cuadrícula

Muestrear de la posterior es una de las técnicas más útiles en la estadística bayesiana. Una vez que tenemos una muestra de valores del parámetro, podemos aproximar cualquier cantidad posterior calculando la función correspondiente en la muestra. Esto convierte integrales difíciles en promedios simples.

3.1.1 ¿Por qué muestrear?

La distribución a posteriori es una distribución de probabilidad sobre los parámetros. En principio, podríamos calcular cualquier cantidad que queramos por integración:

$E[f(\theta)] = \int f(\theta) \, p(\theta \mid D) \, d\theta$

Pero en la práctica, esta integral a menudo es imposible de calcular analíticamente. El muestreo resuelve este problema. Si tenemos $S$ muestras $\theta_1, \theta_2, \ldots, \theta_S$ de la posterior, podemos aproximar la integral por el promedio muestral:

$\hat{E}[f(\theta)] \approx \frac{1}{S} \sum_{s=1}^{S} f(\theta_s)$

Este es el principio de Monte Carlo: reemplazar integración con promedios. La ley de los grandes números garantiza que la aproximación mejora a medida que $S$ aumenta.

3.1.2 Muestreo desde una posterior aproximada por malla

En el Capítulo 2, calculamos la posterior en una malla de valores del parámetro. Cada punto de la malla tiene una probabilidad a posteriori. Para muestrear de esta distribución, simplemente extraemos puntos de la malla con probabilidades proporcionales a sus pesos a posteriori.

# Sample from a grid-approximate posterior
# Globe tossing: 6 waters in 9 tosses

# 1. Build the grid and posterior

p*grid <- seq(from = 0, to = 1, length.out = 1000)
prior <- rep(1, 1000)
likelihood <- dbinom(6, size = 9, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood _ prior)

# 2. Sample 10,000 values from the grid

set.seed(42)
samples <- sample(p_grid, size = 1e4, replace = TRUE, prob = posterior)

# 3. Plot the samples

par(mfrow = c(1, 2))

# Grid posterior

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "p", ylab = "Posterior density",
main = "Grid posterior")

# Sampled posterior

hist(samples, breaks = 50, col = "steelblue", border = "white",
freq = FALSE, xlab = "p",
main = "Sampled posterior (10,000 draws)")
lines(p_grid, posterior * 1000, lwd = 2, col = "firebrick")

par(mfrow = c(1, 1))

cat("Number of samples:", length(samples), "\n")
cat("Sample mean:", round(mean(samples), 4), "\n")
cat("Sample SD:", round(sd(samples), 4), "\n")

Click Run to execute this code in your browser with WebR

La función sample extrae 10.000 valores de la malla, con probabilidades proporcionales a la posterior. Cada extracción es un valor plausible de $p$ . El histograma de estas extracciones aproxima la distribución a posteriori.

3.1.3 ¿Qué podemos hacer con las muestras?

Una vez que tenemos muestras de la posterior, podemos:

Estimar la media a posteriori: mean(samples)
Estimar intervalos de credibilidad: quantile(samples, c(0.05, 0.95))
Estimar probabilidades a posteriori: mean(samples > 0.5)
Simular predicciones: rbinom(1e4, size = 9, prob = samples)
Calcular cualquier función de los parámetros: mean(samples^2) para $E[p^2]$

La idea clave es que toda pregunta sobre la posterior se convierte en una pregunta empírica sobre las muestras. Nunca necesitamos calcular otra integral — solo calculamos sobre las muestras y tomamos el promedio.

# Everything we can do with posterior samples
# Using the samples from the previous code block

# Posterior mean

cat("Posterior mean of p:", round(mean(samples), 4), "\n")

# Posterior median

cat("Posterior median:", round(median(samples), 4), "\n")

# 90% credible interval (percentile)

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

# Posterior probability that p > 0.5

cat("P(p > 0.5):", round(mean(samples > 0.5), 4), "\n")

# Posterior probability that p > 0.6

cat("P(p > 0.6):", round(mean(samples > 0.6), 4), "\n")

# Expected value of p^2

cat("E[p^2]:", round(mean(samples^2), 4), "\n")

# Variance of p

cat("Var(p):", round(var(samples), 6), "\n")

Click Run to execute this code in your browser with WebR

Rethinking

El muestreo es el disolvente universal. La belleza del muestreo es que unifica todos los cálculos a posteriori. Ya sea que quiera una media, un intervalo, una probabilidad o una predicción, la respuesta es la misma: calcule la función sobre las muestras y tome el promedio. Por esto MCMC (Capítulo 8) es tan poderoso — produce muestras de cualquier posterior, no importa cuán compleja. Una vez que tiene muestras, el resto es aritmética.

Overthinking

¿Por qué funciona el muestreo? La ley de los grandes números garantiza que si $\theta_1, \ldots, \theta_S$ son muestras i.i.d. de $p(\theta \mid D)$ , entonces:

$\frac{1}{S} \sum_{s=1}^S f(\theta_s) \xrightarrow{a.s.} E[f(\theta)] = \int f(\theta) \, p(\theta \mid D) \, d\theta$

La tasa de convergencia es $O(1/\sqrt{S})$ , por el teorema del límite central. Esto significa que para reducir a la mitad el error de Monte Carlo, necesitamos cuatro veces más muestras. Con 10.000 muestras, el error es típicamente alrededor del 1% de la desviación estándar — más que suficiente para propósitos prácticos.

Es importante destacar que la tasa de convergencia no depende de la dimensión de $\theta$ . Por esto el muestreo escala a modelos de alta dimensión mientras que la aproximación por malla no. La aproximación por malla requiere $g^k$ puntos para $k$ parámetros, pero el muestreo requiere el mismo número de muestras independientemente de $k$ (aunque la mezcla de la cadena de Markov puede ser más lenta en altas dimensiones).

Resumen de la Sección 3.1

El muestreo reemplaza integración con promedios: $\hat{E}[f(\theta)] \approx \frac{1}{S} \sum f(\theta_s)$ .
Para muestrear de una posterior por malla, extraiga puntos de la malla con probabilidad proporcional a su peso a posteriori.
Una vez que tenemos muestras, toda pregunta posterior se convierte en una pregunta empírica.
La ley de los grandes números garantiza convergencia; la tasa es $O(1/\sqrt{S})$ .
El muestreo escala a altas dimensiones, a diferencia de la aproximación por malla.

Práctica

Muestreo para resumir