4.3 Un modelo gaussiano de altura

Un modelo Gaussiano simple de estatura asume que las estaturas observadas están distribuidas normalmente alrededor de una media con una desviación estándar fija. Usando los datos de Howell, esta sección muestra cómo ajustar tal modelo e interpretar la posterior para la media y la desviación estándar.

4.3.1 Los datos

Los datos de Howell contienen mediciones de estatura y peso del pueblo !Kung San del desierto de Kalahari. Usaremos este conjunto de datos a lo largo del capítulo.

# Simulate the Howell height data (since we cannot load the real dataset in WebR)
# The real data has 544 individuals, adults have mean ~154 cm, SD ~7.7 cm
set.seed(42)

# Simulate adult heights (age >= 18)

n <- 352
height <- rnorm(n, mean = 154.6, sd = 7.7)

# Summary statistics

cat("Howell height data (simulated adults):\n")
cat("n:", length(height), "\n")
cat("Mean:", round(mean(height), 2), "cm\n")
cat("SD:", round(sd(height), 2), "cm\n")
cat("Min:", round(min(height), 2), "cm\n")
cat("Max:", round(max(height), 2), "cm\n")

# Plot

hist(height, breaks = 30, col = "steelblue", border = "white",
freq = FALSE, xlab = "Height (cm)",
main = "Howell height data (adults)")
curve(dnorm(x, mean = mean(height), sd = sd(height)),
add = TRUE, lwd = 2, col = "firebrick")

Click Run to execute this code in your browser with WebR

4.3.2 El modelo

El modelo Gaussiano de estatura es:

$h_i \sim \text{Normal}(\mu, \sigma)$ $\mu \sim \text{Normal}(178, 20)$ $\sigma \sim \text{Uniform}(0, 50)$

La distribución a priori para $\mu$ está centrada en 178 cm (una estatura promedio razonable) con una desviación estándar de 20 cm. La distribución a priori para $\sigma$ es uniforme en $[0, 50]$ — sabemos que la desviación estándar debe ser positiva y menor que 50 cm.

4.3.3 Simulación predictiva a priori

Antes de ajustar el modelo, debemos verificar qué implican las distribuciones a priori sobre los datos. La simulación predictiva a priori extrae de la distribución a priori y genera datos sintéticos. Si los datos sintéticos parecen irrazonables, las distribuciones a priori son demasiado anchas.

# Prior predictive simulation for the height model
set.seed(42)

# Draw from the prior for mu and sigma

n_sim <- 100
mu_prior <- rnorm(n_sim, mean = 178, sd = 20)
sigma_prior <- runif(n_sim, 0, 50)

# Plot the prior predictive distribution

par(mfrow = c(1, 2))

# Prior for mu

hist(mu_prior, breaks = 20, col = "steelblue", border = "white",
main = "Prior for mu", xlab = "mu (cm)")

# Prior predictive: simulate heights from each prior draw

plot(NULL, xlim = c(50, 300), ylim = c(0, 0.08),
xlab = "Height (cm)", ylab = "Density",
main = "Prior predictive: 20 prior draws")
for (i in 1:20) {
curve(dnorm(x, mean = mu_prior[i], sd = sigma_prior[i]),
add = TRUE, col = rgb(0.27, 0.51, 0.71, 0.3))
}
par(mfrow = c(1, 1))

cat("Prior for mu: Normal(178, 20)\n")
cat(" Range:", round(range(mu_prior), 1), "\n")
cat("Prior for sigma: Uniform(0, 50)\n")
cat(" Range:", round(range(sigma_prior), 1), "\n")
cat("\nSome prior predictions are unreasonable (e.g., heights > 250 cm)\n")
cat("But the data will dominate the posterior\n")

Click Run to execute this code in your browser with WebR

La distribución predictiva a priori incluye algunas estaturas irrazonables (e.g., 250 cm), pero la distribución a priori es lo suficientemente amplia para dejar que los datos dominen. Si quisiéramos distribuciones a priori más informativas, podríamos restringirlas (e.g., $\mu \sim \text{Normal}(178, 5)$ ).

4.3.4 Aproximación por malla

Podemos ajustar el modelo usando aproximación por malla. Con dos parámetros ( $\mu$ y $\sigma$ ), la malla es bidimensional.

# Grid approximation for the Gaussian height model
# Using simulated Howell data

set.seed(42)
height <- rnorm(352, mean = 154.6, sd = 7.7)

# Grid for mu and sigma

mu_grid <- seq(140, 170, length.out = 100)
sigma_grid <- seq(4, 12, length.out = 100)

# Compute log-likelihood on the grid

# For efficiency, compute on a coarse grid

log_lik <- matrix(0, nrow = length(mu_grid), ncol = length(sigma_grid))
for (i in 1:length(mu_grid)) {
for (j in 1:length(sigma_grid)) {
log_lik[i, j] <- sum(dnorm(height, mean = mu_grid[i], sd = sigma_grid[j], log = TRUE))
}
}

# Add log-prior

log_prior_mu <- dnorm(mu_grid, 178, 20, log = TRUE)
log_prior_sigma <- dunif(sigma_grid, 0, 50, log = TRUE)
log_post <- log_lik + outer(log_prior_mu, log_prior_sigma)

# Normalize

post <- exp(log_post - max(log_post))
post <- post / sum(post)

# Marginal posteriors

post_mu <- rowSums(post)
post_sigma <- colSums(post)

# Plot

par(mfrow = c(1, 2))
plot(mu_grid, post_mu, type = "l", lwd = 2, col = "steelblue",
xlab = "mu", ylab = "Posterior density", main = "Posterior for mu")
abline(v = mean(height), lty = 2, col = "firebrick")
plot(sigma_grid, post_sigma, type = "l", lwd = 2, col = "firebrick",
xlab = "sigma", ylab = "Posterior density", main = "Posterior for sigma")
abline(v = sd(height), lty = 2, col = "steelblue")
par(mfrow = c(1, 1))

cat("Posterior mean of mu:", round(sum(mu*grid * post*mu), 2), "\n")
cat("Posterior mean of sigma:", round(sum(sigma_grid * post_sigma), 2), "\n")
cat("Sample mean:", round(mean(height), 2), "\n")
cat("Sample SD:", round(sd(height), 2), "\n")

Click Run to execute this code in your browser with WebR

4.3.5 Aproximación cuadrática

Para este modelo, la posterior es aproximadamente Gaussiana, así que la aproximación cuadrática funciona bien.

# Quadratic approximation for the Gaussian height model
# The MAP estimate is the sample mean and sample SD

set.seed(42)
height <- rnorm(352, mean = 154.6, sd = 7.7)

# MAP estimates

mu_map <- mean(height)
sigma_map <- sd(height)

# Standard errors (from the Hessian)

# For the normal model: SE(mu) = sigma / sqrt(n), SE(sigma) = sigma / sqrt(2n)

n <- length(height)
se_mu <- sigma_map / sqrt(n)
se_sigma <- sigma_map / sqrt(2 * n)

cat("Quadratic approximation:\n")
cat("mu: ", round(mu_map, 2), " +/- ", round(se_mu, 3), "\n")
cat("sigma: ", round(sigma_map, 2), " +/- ", round(se_sigma, 3), "\n")
cat("\n95% interval for mu: [", round(mu_map - 1.96 * se_mu, 2), ",",
round(mu_map + 1.96 * se_mu, 2), "]\n")
cat("95% interval for sigma: [", round(sigma_map - 1.96 * se_sigma, 2), ",",
round(sigma_map + 1.96 * se_sigma, 2), "]\n")

Click Run to execute this code in your browser with WebR

Rethinking

Siempre verifique sus distribuciones a priori. La simulación predictiva a priori es el hábito más importante en el modelado bayesiano. Antes de ajustar el modelo, simule desde las distribuciones a priori y verifique si los datos implícitos son razonables. Si su distribución a priori dice que estaturas de 300 cm son plausibles, quizás quiera restringirla. Si su distribución a priori dice que las estaturas deben estar entre 150 y 160 cm, quizás quiera aflojarla. El objetivo es una distribución a priori débilmente informativa: que descarte valores absurdos pero deje que los datos hablen.

Overthinking

La distribución a priori de Jeffreys para sigma. La distribución uniforme sobre $\sigma$ no es la elección más natural. La distribución a priori de Jeffreys para la desviación estándar de una distribución normal es $p(\sigma) \propto 1/\sigma$ , que es invariante a escala: no depende de las unidades de medición. Usando la distribución a priori de Jeffreys, la posterior para $\sigma$ es:

$\sigma \mid D \sim \text{Inverse-Gamma}\left(\frac{n-1}{2}, \frac{(n-1)s^2}{2}\right)$

donde $s^2$ es la varianza muestral. Para $n$ grande, esto es muy cercano al resultado con la distribución uniforme. La diferencia importa más para muestras pequeñas.

Resumen de la Sección 4.3

El modelo Gaussiano de estatura tiene dos parámetros: $\mu$ (media) y $\sigma$ (desviación estándar).
La simulación predictiva a priori verifica qué implican las distribuciones a priori sobre los datos antes de ajustar.
La aproximación por malla funciona para 2 parámetros pero es computacionalmente costosa.
La aproximación cuadrática es rápida y precisa para este modelo porque la posterior es casi Gaussiana.
Siempre use distribuciones a priori débilmente informativas que descarten valores absurdos pero dejen que los datos hablen.

Un lenguaje para describir modelos

Añadiendo un predictor