1. Purpose of the Example
The passage presents a hierarchical normal model to illustrate Gibbs sampling and Metropolis algorithm techniques in Bayesian computation.
The goal is not to do a complete applied data analysis but to show how simulation-based computation works in a relatively simple model.
2. The Data
A small experimental dataset contains observations grouped into J groups, where group j has $n_j$ data points:
- Data within each group:
- $y_{ij},\ i = 1,\dots, n_j,\ j = 1,\dots, J$
Total observations:
$n = \sum_{j=1}^J n_j.$
3. The Hierarchical Normal Model
Level 1: Data model
Within each group j:
$y_{ij} \sim N(\theta_j,\ \sigma^2).$
All groups share a common variance σ².
Level 2: Group-level prior
Group means $θ_j$ follow a normal distribution:
$\theta_j \sim N(\mu,\ \tau^2)$
Priors for hyperparameters
A uniform prior on:
$(\mu,\ \log \sigma,\ \log \tau), \qquad \sigma > 0,\ \tau > 0$
which implies:
$p(\mu, \log\sigma, \log\tau) \propto \tau.$
A uniform prior on log τ would make the posterior improper, so they avoid that.
4. Joint Posterior Distribution
The full joint posterior is:
$p(\theta,\ \mu,\ \log\sigma,\ \log\tau\mid y) \ \propto\ \tau \prod_{j=1}^J N(\theta_j \mid \mu,\ \tau^2) \prod_{j=1}^J \ \prod_{i=1}^{n_j} N(y_{ij} \mid \theta_j,\ \sigma^2)$
This is the density that the Gibbs sampler and Metropolis algorithm aim to explore.
5. Overdispersed Starting Points
To initialize the chains:
- For each group j: choose $θ_j$ by drawing a random observation from that group.
- Set μ as the average of the initial $θ_j$’s.
- σ and τ do not need initial values—they are drawn in the first Gibbs steps.
This ensures chains begin from “overdispersed” starting points, improving convergence diagnostics.
6. Gibbs Sampler
All full conditional distributions have conjugate forms.
6.1 Conditional Posterior of $θ_j$
Given μ, σ, τ, y:
$\theta_j \mid \mu,\sigma,\tau,y \sim N(\hat{\theta}_j, V_{\theta_j})$
where
$\hat{\theta}_j = \frac{\frac{1}{\tau^2}\mu + \frac{n_j}{\sigma^2} \bar{y}_j} {\frac{1}{\tau^2} + \frac{n_j}{\sigma^2}}$
$V_{\theta_j} = \frac{1}{\frac{1}{\tau^2} + \frac{n_j}{\sigma^2}}$
The $θ_j$’s are independent conditional on hyperparameters.
6.2 Conditional Posterior of μ
Given θ:
$\mu \mid \theta,\sigma,\tau,y \sim N(\hat{\mu},\ \tau^2/J)$
with
$\hat{\mu} = \frac{1}{J}\sum_{j=1}^J \theta_j$
6.3 Conditional Posterior of σ²
There are n observations, each with mean θ_j. Thus:
$\sigma^2 \mid \theta,\mu,\tau,y \sim \text{Inverse-}\chi^2\left( n,\ \hat{\sigma}^2 \right)$
where
$\hat{\sigma}^2 = \frac{1}{n}\sum_{j=1}^J \sum_{i=1}^{n_j} (y_{ij} – \theta_j)^2$
6.4 Conditional Posterior of τ²
$\tau^2 \mid \theta,\mu,\sigma,y \sim \text{Inverse-}\chi^2\left( J-1,\ \hat{\tau}^2 \right)$
where
$\hat{\tau}^2 = \frac{1}{J – 1}\sum_{j=1}^J (\theta_j – \mu)^2$
The degrees of freedom are $J – 1$ because prior $p(\tau)\propto 1$, not proportional to τ⁻¹.
7. Numerical Results for the Coagulation Data
- They ran 10 parallel Gibbs chains.
- Each chain: 100 iterations.
- Convergence diagnostic (potential scale reduction $\hat{R}$ was near 1.0 for all parameters.
- Results summarized in the below.
Typical findings:
- σ (within-group SD) is estimated reasonably well.
- τ (between-group variation) has higher uncertainty because there are few groups, which is typical.
8. Metropolis Algorithm Approach
Instead of sampling all parameters jointly, they:
- Jump only in the 3-dimensional space: $(\mu,\ \log\sigma,\ \log\tau)$
- Then draw θ from its normal conditional posterior.
Efficient Proposal Distribution
They use a multivariate normal proposal, centered at current point, with covariance equal to the normal approximation of the posterior multiplied by:
$\frac{2.4^2}{d}$
where $d = 3$.
This is a standard optimal scaling result for Metropolis in multidimensional normal targets.
9. Results for Metropolis Algorithm
- 10 parallel chains, 500 iterations each.
- Convergence: $\hat{R} < 1.1$ for all parameters.
- Posterior results matched those from the Gibbs sampler.
- Acceptance rate: 0.35, which is close to the theoretical optimum for dimension 3 under the 2.4/√d rule.
10. Main Takeaways
- The hierarchical normal model provides clean conjugate forms for Gibbs sampling.
- Group-level parameters ($θ_j$) and hyperparameters (μ, σ, τ) all have straightforward full conditionals.
- Gibbs sampler converges quickly due to conjugacy.
- Metropolis algorithm can efficiently explore only the hyperparameter space, then draw θ from conditional normals.
- The uncertainty in τ (between-group variance) remains high with few groups—typical in hierarchical modeling.
- Both sampling methods produce consistent results when tuned properly.