2.4 Haciendo que el modelo funcione

Una vez que un modelo está especificado — verosimilitud, parámetros y distribución a priori — necesitamos calcular la posterior. Este es el corazón práctico de la estadística bayesiana. En esta sección, exploramos tres estrategias para calcular la posterior: aproximación por malla, aproximación cuadrática y cadenas de Markov Monte Carlo.

2.4.1 Aproximación por malla

La aproximación por malla es el método más simple. Disponemos una malla de posibles valores del parámetro, calculamos la posterior en cada punto de la malla y normalizamos. Esto es exactamente lo que hicimos en la Sección 2.1.

La ventaja de la aproximación por malla es que es exacta (hasta la resolución de la malla) y transparente — podemos ver cada paso. La desventaja es que escala mal: con $k$ parámetros y $g$ puntos de malla por parámetro, necesitamos $g^k$ cálculos. Para un modelo con 10 parámetros y 100 puntos de malla cada uno, eso son $100^{10} = 10^{20}$ cálculos — inviable.

# Grid approximation for the globe tossing problem
# Data: 6 waters in 9 tosses

# 1. Define the grid

p_grid <- seq(from = 0, to = 1, length.out = 100)

# 2. Define the prior (uniform)

prior <- rep(1, 100)

# 3. Compute the likelihood

likelihood <- dbinom(6, size = 9, prob = p_grid)

# 4. Compute the unstandardized posterior

unstd_posterior <- likelihood * prior

# 5. Standardize

posterior <- unstd_posterior / sum(unstd_posterior)

# 6. Plot

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "Probability of water (p)",
ylab = "Posterior probability",
main = "Grid approximation: 100 points")
polygon(c(p_grid, rev(p_grid)),
c(posterior, rep(0, length(posterior))),
col = rgb(0.27, 0.51, 0.71, 0.3), border = NA)

# 7. Summarize

cat("Posterior mean:", round(sum(p_grid * posterior), 4), "\n")
cat("Posterior mode:", round(p_grid[which.max(posterior)], 4), "\n")

Click Run to execute this code in your browser with WebR

La aproximación por malla funciona bien para modelos con 1-2 parámetros. Más allá de eso, necesitamos métodos más eficientes.

2.4.2 Aproximación cuadrática

La aproximación cuadrática (también llamada aproximación de Laplace) aproxima la posterior con una distribución Gaussiana multivariada. La idea es encontrar el pico de la posterior (la estimación máxima a posteriori o MAP) y luego aproximar la posterior como una Gaussiana centrada en el pico, con covarianza determinada por la curvatura del logaritmo de la posterior.

El logaritmo de la posterior cerca del pico es aproximadamente cuadrático:

$\log p(\theta \mid D) \approx \log p(\hat{\theta} \mid D) - \frac{1}{2}(\theta - \hat{\theta})^T H (\theta - \hat{\theta})$

donde $\hat{\theta}$ es la estimación MAP y $H$ es la matriz Hessiana (matriz de segundas derivadas) del logaritmo de la posterior. La posterior aproximada es:

$p(\theta \mid D) \approx \text{Normal}(\hat{\theta}, H^{-1})$

# Quadratic approximation for the globe tossing problem
# Compare with grid approximation

# Grid approximation (reference)

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

# Quadratic approximation

# The MAP estimate for a binomial with uniform prior is w/n

w <- 6; n <- 9
p_map <- w / n

# The standard deviation of the quadratic approximation

# For binomial: sigma = sqrt(p*(1-p)/n)

sigma_quad <- sqrt(p_map * (1 - p_map) / n)

cat("MAP estimate:", round(p_map, 4), "\n")
cat("Quadratic approximation sigma:", round(sigma_quad, 4), "\n")

# Compare grid and quadratic

plot(p_grid, post_grid * 1000, type = "l", lwd = 2, col = "steelblue",
xlab = "Probability of water (p)",
ylab = "Density",
main = "Grid vs Quadratic approximation")
curve(dnorm(x, mean = p_map, sd = sigma_quad), from = 0, to = 1,
add = TRUE, lwd = 2, col = "firebrick", lty = 2)
legend("topright",
legend = c("Grid approximation", "Quadratic approximation"),
col = c("steelblue", "firebrick"), lty = c(1, 2), lwd = 2)

Click Run to execute this code in your browser with WebR

La aproximación cuadrática es muy cercana a la aproximación por malla cuando la posterior tiene forma aproximadamente de campana. La ventaja de la aproximación cuadrática es que es rápida — solo necesitamos encontrar el pico y calcular la curvatura. La desventaja es que falla cuando la posterior es asimétrica, multimodal o tiene soporte acotado (como $p \in [0, 1]$ ).

2.4.3 Cadenas de Markov Monte Carlo

Cadenas de Markov Monte Carlo (MCMC) es el método más flexible. En lugar de calcular la posterior en una malla o aproximarla con una Gaussiana, MCMC extrae muestras de la posterior. Estas muestras pueden usarse para calcular cualquier cantidad de interés (medias, intervalos, predicciones) tomando promedios.

La idea clave es construir una cadena de Markov — una secuencia de valores del parámetro donde cada valor depende solo del anterior — cuya distribución estacionaria es la posterior. Después de que la cadena ha convergido, las muestras están distribuidas según la posterior.

MCMC escala a modelos con cientos o miles de parámetros, haciéndolo el caballo de batalla de la estadística bayesiana moderna. Estudiaremos MCMC en detalle en el Capítulo 8.

# MCMC from scratch: Metropolis algorithm for the globe tossing problem
# Data: 6 waters in 9 tosses

set.seed(42)
n_samples <- 1000
p_current <- 0.5 # starting value
samples <- numeric(n_samples)

for (i in 1:n_samples) {

# Propose a new value by adding a small random jump

p_proposal <- p_current + runif(1, -0.1, 0.1)

# Ensure the proposal stays in [0, 1]

if (p_proposal < 0) p_proposal <- abs(p_proposal)
if (p_proposal > 1) p_proposal <- 2 - p_proposal

# Compute the acceptance ratio

# Posterior is proportional to likelihood * prior

# Prior is uniform, so it cancels

log_post_current <- dbinom(6, 9, p_current, log = TRUE)
log_post_proposal <- dbinom(6, 9, p_proposal, log = TRUE)

# Accept or reject

if (log(runif(1)) < (log_post_proposal - log_post_current)) {
p_current <- p_proposal
}
samples[i] <- p_current
}

# Plot the samples

par(mfrow = c(1, 2))

# Trace plot

plot(1:n_samples, samples, type = "l", col = "steelblue",
xlab = "Iteration", ylab = "p",
main = "MCMC trace plot")

# Histogram of samples

hist(samples, breaks = 30, col = "steelblue", border = "white",
freq = FALSE, xlab = "p",
main = "MCMC posterior samples")

# Overlay grid approximation

p_grid <- seq(0, 1, length.out = 100)
post_grid <- dbinom(6, 9, p_grid) / sum(dbinom(6, 9, p_grid))
lines(p_grid, post_grid * 100, lwd = 2, col = "firebrick")

par(mfrow = c(1, 1))

cat("MCMC posterior mean:", round(mean(samples), 4), "\n")
cat("MCMC posterior SD:", round(sd(samples), 4), "\n")
cat("MCMC 95% interval: [", round(quantile(samples, 0.025), 4), ",",
round(quantile(samples, 0.975), 4), "]\n")

Click Run to execute this code in your browser with WebR

El gráfico de traza muestra la cadena explorando el espacio de parámetros. El histograma de muestras aproxima la distribución a posteriori. Con más muestras, la aproximación mejora. MCMC es el método más poderoso, pero requiere un diagnóstico cuidadoso (Capítulo 8).

2.4.4 Comparando los tres métodos

| Método | ¿Exacto? | ¿Escalable? | ¿Rápido? | Cuándo usar | | ---------- | ------------------- | ----------- | ------------------ | ---------------------------------- | | Malla | Sí (hasta la malla) | No ( $g^k$ ) | Lento para $k > 2$ | Enseñanza, 1-2 parámetros | | Cuadrático | Aproximado | Sí | Muy rápido | Posteriors casi Gaussianas | | MCMC | Asintóticamente | Sí | Moderado | Modelos complejos, cualquier forma |

Rethinking

El mundo pequeño y el mundo grande. Los tres métodos operan en el mundo pequeño — el modelo matemático que hemos construido. El mundo pequeño es autoconsistente: dado el modelo, la posterior está bien definida, y los métodos la aproximan. Pero el mundo grande — la realidad — puede no coincidir con el modelo. Un cálculo perfecto de la posterior equivocada sigue siendo equivocado. Los métodos de esta sección nos dicen cómo calcular la posterior con precisión. No nos dicen si el modelo es correcto. Ese es el trabajo de la crítica del modelo (Capítulo 6) y del juicio científico.

Overthinking

La Hessiana y la aproximación cuadrática. La matriz Hessiana $H$ del logaritmo de la posterior contiene las segundas derivadas:

$H_{ij} = \frac{\partial^2 \log p(\theta \mid D)}{\partial \theta_i \partial \theta_j}$

La inversa de la Hessiana $H^{-1}$ es la matriz de covarianza de la aproximación cuadrática. Para un modelo con un solo parámetro, $H$ es un escalar y $\sigma = 1/\sqrt{-H}$ . Para el modelo binomial con $w$ éxitos en $n$ ensayos y una distribución uniforme a priori:

$\log p(p \mid D) = w \log p + (n-w) \log(1-p) + \text{const}$

$\frac{d^2 \log p}{dp^2} = -\frac{w}{p^2} - \frac{n-w}{(1-p)^2}$

En el MAP $\hat{p} = w/n$ :

$H = -\frac{w}{(w/n)^2} - \frac{n-w}{(1-w/n)^2} = -\frac{n^2}{w} - \frac{n^2}{n-w} = -\frac{n^3}{w(n-w)}$

$\sigma = \frac{1}{\sqrt{n^3 / (w(n-w))}} = \sqrt{\frac{w(n-w)}{n^3}} = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$

Este es el familiar error estándar de una proporción, confirmando que la aproximación cuadrática recupera el resultado clásico.

Resumen de la Sección 2.4

La aproximación por malla calcula la posterior en valores discretos del parámetro. Exacta pero escala como $g^k$ .
La aproximación cuadrática ajusta una Gaussiana al pico de la posterior. Rápida y precisa para posteriors con forma de campana.
MCMC extrae muestras de la posterior. Flexible y escala a modelos complejos, pero requiere diagnóstico.
Los tres métodos intercambian precisión, velocidad y escalabilidad.
Todos los métodos operan en el mundo pequeño del modelo — no garantizan que el modelo sea correcto.

Componentes del modelo

Resumen