4.3 A Gaussian model of height

A simple Gaussian model of height assumes that observed heights are normally distributed around a mean with a fixed standard deviation. Using the Howell data, this section shows how to fit such a model and interpret the posterior for the mean and standard deviation.

4.3.1 The data

The Howell data contains height and weight measurements from the !Kung San people of the Kalahari Desert. We will use this dataset throughout the chapter.

# 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 The model

The Gaussian model of height is:

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

The prior for $\mu$ is centered at 178 cm (a reasonable average height) with a standard deviation of 20 cm. The prior for $\sigma$ is uniform on $[0, 50]$ — we know the standard deviation must be positive and less than 50 cm.

4.3.3 Prior predictive simulation

Before fitting the model, we should check what the priors imply about the data. Prior predictive simulation draws from the prior and generates synthetic data. If the synthetic data looks unreasonable, the priors are too wide.

# 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

The prior predictive distribution includes some unreasonable heights (e.g., 250 cm), but the prior is wide enough to let the data dominate. If we wanted more informative priors, we could tighten them (e.g., $\mu \sim \text{Normal}(178, 5)$ ).

4.3.4 Grid approximation

We can fit the model using grid approximation. With two parameters ( $\mu$ and $\sigma$ ), the grid is 2-dimensional.

# 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 Quadratic approximation

For this model, the posterior is approximately Gaussian, so quadratic approximation works well.

# 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

Always check your priors. Prior predictive simulation is the most important habit in Bayesian modeling. Before fitting the model, simulate from the priors and check whether the implied data are reasonable. If your prior says that heights of 300 cm are plausible, you might want to tighten it. If your prior says that heights must be between 150 and 160 cm, you might want to loosen it. The goal is a prior that is weakly informative: it rules out absurd values but lets the data speak.

Overthinking

The Jeffreys prior for sigma. The uniform prior on $\sigma$ is not the most natural choice. The Jeffreys prior for the standard deviation of a normal distribution is $p(\sigma) \propto 1/\sigma$ , which is scale-invariant: it does not depend on the units of measurement. Using the Jeffreys prior, the posterior for $\sigma$ is:

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

where $s^2$ is the sample variance. For large $n$ , this is very close to the result with the uniform prior. The difference matters most for small samples.

Summary of Section 4.3

The Gaussian model of height has two parameters: $\mu$ (mean) and $\sigma$ (standard deviation).
Prior predictive simulation checks what the priors imply about the data before fitting.
Grid approximation works for 2 parameters but is computationally expensive.
Quadratic approximation is fast and accurate for this model because the posterior is nearly Gaussian.
Always use weakly informative priors that rule out absurd values but let the data speak.

A language for describing models

Adding a predictor