Complete steps of the Bayesian Hierarchical Normal Model

Introduction

This document provides a complete, step-by-step algebraic derivation of the Bayesian hierarchical model. We start from the joint posterior distribution and show every algebraic step to reach the final factorization.

The Model

The hierarchical model is specified as:

\[ \begin{aligned} y_i &\mid \mu, \sigma^2 \sim N(\mu, \sigma^2), \quad i = 1, 2, \ldots, n \\ \mu &\mid \sigma^2 \sim N\left(\theta, \frac{\sigma^2}{n_0}\right) \\ \sigma^2 &\sim IG(a_0, b_0) \end{aligned} \]

where \(\theta\), \(n_0\), \(a_0\), and \(b_0\) are known scalars.

---

Part 1: From Bayes’ Theorem to the Joint Posterior

Step 1: Apply Bayes’ Theorem

By Bayes’ theorem, the joint posterior is proportional to the prior times the likelihood:

\[ p(\mu, \sigma^2 \mid y) \propto p(\mu, \sigma^2) \times p(y \mid \mu, \sigma^2) \]

Step 2: Factor the Joint Prior

Using the hierarchical structure:

\[ p(\mu, \sigma^2) = p(\mu \mid \sigma^2) \, p(\sigma^2) \]

Step 3: Write Out All Components

\[ p(\mu, \sigma^2 \mid y) \propto \underbrace{IG(\sigma^2 \mid a_0, b_0)}_{\text{prior for } \sigma^2} \times \underbrace{N\left(\mu \mid \theta, \frac{\sigma^2}{n_0}\right)}_{\text{conditional prior for } \mu} \times \underbrace{\prod_{i=1}^n N(y_i \mid \mu, \sigma^2)}_{\text{likelihood}} \tag{A1} \]

---

Part 2: Simplifying the Densities

Step 4: Write the Inverse-Gamma Density

The Inverse-Gamma density with shape \(a_0\) and scale \(b_0\) is:

\[ IG(\sigma^2 \mid a_0, b_0) = \frac{b_0^{a_0}}{\Gamma(a_0)} (\sigma^2)^{-(a_0+1)} \exp\left(-\frac{b_0}{\sigma^2}\right) \]

Key point: When we use \(\propto\), we drop any terms that don’t depend on \(\mu\) or \(\sigma^2\). So:

\[ IG(\sigma^2 \mid a_0, b_0) \propto (\sigma^2)^{-(a_0+1)} \exp\left(-\frac{b_0}{\sigma^2}\right) \tag{A2} \]

Step 5: Write the Conditional Prior for \(\mu\)

The Normal density for \(\mu\) given \(\sigma^2\):

\[ N\left(\mu \mid \theta, \frac{\sigma^2}{n_0}\right) = \sqrt{\frac{n_0}{2\pi\sigma^2}} \exp\left(-\frac{n_0(\mu-\theta)^2}{2\sigma^2}\right) \]

Dropping constants:

\[ N\left(\mu \mid \theta, \frac{\sigma^2}{n_0}\right) \propto (\sigma^2)^{-1/2} \exp\left(-\frac{n_0(\mu-\theta)^2}{2\sigma^2}\right) \tag{A3} \]

Step 6: Write the Likelihood

Each observation has density:

\[ N(y_i \mid \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(y_i-\mu)^2}{2\sigma^2}\right) \]

The product over \(i=1,\ldots,n\):

\[ \prod_{i=1}^n N(y_i \mid \mu, \sigma^2) = (2\pi\sigma^2)^{-n/2} \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i-\mu)^2\right) \]

Dropping constants:

\[ \prod_{i=1}^n N(y_i \mid \mu, \sigma^2) \propto (\sigma^2)^{-n/2} \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i-\mu)^2\right) \tag{A4} \]

---

Part 3: Combining All Components

Step 7: Multiply Everything Together

Substitute (A2), (A3), and (A4) into (A1):

\[ p(\mu, \sigma^2 \mid y) \propto \left[(\sigma^2)^{-(a_0+1)} \exp\left(-\frac{b_0}{\sigma^2}\right)\right] \times \left[(\sigma^2)^{-1/2} \exp\left(-\frac{n_0(\mu-\theta)^2}{2\sigma^2}\right)\right] \times \left[(\sigma^2)^{-n/2} \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i-\mu)^2\right)\right] \tag{A5} \]

Step 8: Combine the Powers of \(\sigma^2\)

When multiplying terms with the same base, we add the exponents:

\[ (\sigma^2)^{-(a_0+1)} \times (\sigma^2)^{-1/2} \times (\sigma^2)^{-n/2} = (\sigma^2)^{-\left[(a_0+1) + \frac{1}{2} + \frac{n}{2}\right]} \]

Simplify the exponent:

\[ (a_0+1) + \frac{1}{2} + \frac{n}{2} = a_0 + 1 + \frac{n+1}{2} \]

Therefore:

\[ (\sigma^2)^{-\left(a_0 + 1 + \frac{n+1}{2}\right)} \tag{A6} \]

Step 9: Combine the Exponentials

When multiplying exponentials, we add the exponents:

\[ \exp\left(-\frac{b_0}{\sigma^2}\right) \times \exp\left(-\frac{n_0(\mu-\theta)^2}{2\sigma^2}\right) \times \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i-\mu)^2\right) \]

This equals:

\[ \exp\left[ -\frac{b_0}{\sigma^2} - \frac{n_0(\mu-\theta)^2}{2\sigma^2} - \frac{1}{2\sigma^2} \sum_{i=1}^n (y_i-\mu)^2 \right] \]

Step 10: Factor Out \(-\frac{1}{\sigma^2}\)

Factor \(-\frac{1}{\sigma^2}\) from all terms in the exponent:

\[ \exp\left[ -\frac{1}{\sigma^2} \left( b_0 + \frac{n_0(\mu-\theta)^2}{2} + \frac{1}{2} \sum_{i=1}^n (y_i-\mu)^2 \right) \right] \tag{A7} \]

Step 11: Define \(Q(\mu, \theta, y)\)

Define:

\[ Q(\mu, \theta, y) = n_0(\mu-\theta)^2 + \sum_{i=1}^n (y_i-\mu)^2 \tag{A8} \]

Then:

\[ \frac{n_0(\mu-\theta)^2}{2} + \frac{1}{2} \sum_{i=1}^n (y_i-\mu)^2 = \frac{1}{2} Q(\mu, \theta, y) \]

Step 12: Final Joint Posterior Form

Combining (A6) and (A7) with the definition of \(Q\):

\[ \boxed{ p(\mu, \sigma^2 \mid y) \propto (\sigma^2)^{-\left(a_0 + 1 + \frac{n+1}{2}\right)} \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{1}{2}Q(\mu, \theta, y)\right)\right] } \tag{A9} \]

---

Part 4: Completing the Square in \(\mu\)

Step 13: Expand \(Q(\mu, \theta, y)\)

Recall:

\[ Q(\mu, \theta, y) = n_0(\mu-\theta)^2 + \sum_{i=1}^n (y_i-\mu)^2 \tag{A10} \]

Step 13.1: Expand \(\sum_{i=1}^n (y_i-\mu)^2\)

First, rewrite the sum using the sample mean \(\bar{y}\):

\[ \sum_{i=1}^n (y_i-\mu)^2 = \sum_{i=1}^n [(y_i - \bar{y}) + (\bar{y} - \mu)]^2 \]

Expand:

\[ = \sum_{i=1}^n (y_i - \bar{y})^2 + 2(\bar{y} - \mu)\sum_{i=1}^n (y_i - \bar{y}) + n(\bar{y} - \mu)^2 \]

But \(\sum_{i=1}^n (y_i - \bar{y}) = 0\), so:

\[ \sum_{i=1}^n (y_i-\mu)^2 = \sum_{i=1}^n (y_i - \bar{y})^2 + n(\bar{y} - \mu)^2 \]

And \(\sum_{i=1}^n (y_i - \bar{y})^2 = (n-1)s^2\), where \(s^2\) is the sample variance.

So:

\[ \sum_{i=1}^n (y_i-\mu)^2 = (n-1)s^2 + n(\bar{y} - \mu)^2 \tag{A11} \]

Step 13.2: Write \(Q(\mu, \theta, y)\)

Substitute (A11) into (A10):

\[ Q(\mu, \theta, y) = n_0(\mu-\theta)^2 + (n-1)s^2 + n(\bar{y} - \mu)^2 \tag{A12} \]

Step 14: Complete the Square in \(\mu\)

We want to collect all \(\mu\) terms into a single quadratic. Expand the two quadratic terms:

\[ n_0(\mu-\theta)^2 + n(\bar{y} - \mu)^2 \]

Step 14.1: Expand each term

\[ n_0(\mu-\theta)^2 = n_0(\mu^2 - 2\mu\theta + \theta^2) = n_0\mu^2 - 2n_0\theta\mu + n_0\theta^2 \]

\[ n(\bar{y} - \mu)^2 = n(\mu^2 - 2\mu\bar{y} + \bar{y}^2) = n\mu^2 - 2n\bar{y}\mu + n\bar{y}^2 \]

Step 14.2: Add them together

\[ n_0(\mu-\theta)^2 + n(\bar{y} - \mu)^2 = (n_0 + n)\mu^2 - 2(n_0\theta + n\bar{y})\mu + (n_0\theta^2 + n\bar{y}^2) \tag{A13} \]

Step 14.3: Complete the square

We want to write (A13) in the form:

\[ (n_0 + n)\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2 + \text{constant} \]

Expand the proposed form:

\[ (n_0 + n)\left(\mu^2 - 2\mu\frac{n_0\theta + n\bar{y}}{n_0 + n} + \left(\frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2\right) \]

\[ = (n_0 + n)\mu^2 - 2(n_0\theta + n\bar{y})\mu + \frac{(n_0\theta + n\bar{y})^2}{n_0 + n} \]

This matches (A13) except the constant term. So:

\[ n_0(\mu-\theta)^2 + n(\bar{y} - \mu)^2 = (n_0 + n)\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2 + C \]

where:

\[ C = n_0\theta^2 + n\bar{y}^2 - \frac{(n_0\theta + n\bar{y})^2}{n_0 + n} \]

Step 14.4: Simplify the constant \(C\)

Let’s simplify \(C\):

\[ C = n_0\theta^2 + n\bar{y}^2 - \frac{n_0^2\theta^2 + 2n_0n\theta\bar{y} + n^2\bar{y}^2}{n_0 + n} \]

Put the first two terms over common denominator \(n_0 + n\):

\[ C = \frac{(n_0\theta^2 + n\bar{y}^2)(n_0 + n) - (n_0^2\theta^2 + 2n_0n\theta\bar{y} + n^2\bar{y}^2)}{n_0 + n} \]

Expand the numerator:

\[ = \frac{n_0^2\theta^2 + n_0n\theta^2 + n_0n\bar{y}^2 + n^2\bar{y}^2 - n_0^2\theta^2 - 2n_0n\theta\bar{y} - n^2\bar{y}^2}{n_0 + n} \]

Cancel terms:

\[ = \frac{n_0n\theta^2 + n_0n\bar{y}^2 - 2n_0n\theta\bar{y}}{n_0 + n} \]

\[ = \frac{n_0n(\theta^2 - 2\theta\bar{y} + \bar{y}^2)}{n_0 + n} \]

\[ = \frac{n_0 n}{n_0 + n}(\bar{y} - \theta)^2 \tag{A14} \]

Step 14.5: Final completed square form

So (A13) becomes:

\[ n_0(\mu-\theta)^2 + n(\bar{y} - \mu)^2 = (n_0 + n)\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2 + \frac{n_0 n}{n_0 + n}(\bar{y} - \theta)^2 \tag{A15} \]

Step 15: Rewrite \(Q(\mu, \theta, y)\)

Substitute (A15) into (A12):

\[ Q(\mu, \theta, y) = (n_0 + n)\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2 + \frac{n_0 n}{n_0 + n}(\bar{y} - \theta)^2 + (n-1)s^2 \]

Define:

\[ f(\theta, \bar{y}) = (n-1)s^2 + \frac{n_0 n}{n_0 + n}(\bar{y} - \theta)^2 \tag{A16} \]

Then:

\[ Q(\mu, \theta, y) = (n_0 + n)\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2 + f(\theta, \bar{y}) \tag{A17} \]

---

Part 5: Factorizing the Joint Posterior

Step 16: Substitute \(Q\) into the Joint Posterior

Starting from (A9):

\[ p(\mu, \sigma^2 \mid y) \propto (\sigma^2)^{-\left(a_0 + 1 + \frac{n+1}{2}\right)} \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{1}{2}Q(\mu, \theta, y)\right)\right] \]

Substitute (A17):

\[ p(\mu, \sigma^2 \mid y) \propto (\sigma^2)^{-\left(a_0 + 1 + \frac{n+1}{2}\right)} \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{1}{2}\left((n_0 + n)\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2 + f(\theta, \bar{y})\right)\right)\right] \tag{A18} \]

Step 17: Split the Exponential

Separate the terms inside the exponential:

\[ \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{f(\theta, \bar{y})}{2} + \frac{n_0 + n}{2}\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2\right)\right] \]

\[ = \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{f(\theta, \bar{y})}{2}\right)\right] \times \exp\left[-\frac{n_0 + n}{2\sigma^2}\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2\right] \tag{A19} \]

Step 18: Split the Power of \(\sigma^2\)

We need to split the power of \(\sigma^2\) into two parts that match the Inverse-Gamma and Normal forms:

\[ (\sigma^2)^{-\left(a_0 + 1 + \frac{n+1}{2}\right)} = (\sigma^2)^{-\left(a_0 + \frac{n}{2} + 1\right)} \times (\sigma^2)^{-1/2} \]

Check: \(a_0 + \frac{n}{2} + 1 + \frac{1}{2} = a_0 + 1 + \frac{n+1}{2}\) ✓

So:

\[ (\sigma^2)^{-\left(a_0 + 1 + \frac{n+1}{2}\right)} = (\sigma^2)^{-\left(a_0 + \frac{n}{2} + 1\right)} \times (\sigma^2)^{-1/2} \tag{A20} \]

Step 19: Put Everything Together

Substitute (A19) and (A20) into (A18):

\[ p(\mu, \sigma^2 \mid y) \propto \underbrace{ (\sigma^2)^{-\left(a_0 + \frac{n}{2} + 1\right)} \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{f(\theta, \bar{y})}{2}\right)\right] }_{\text{Inverse-Gamma part}} \times \underbrace{ (\sigma^2)^{-1/2} \exp\left[-\frac{n_0 + n}{2\sigma^2}\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2\right] }_{\text{Normal part}} \tag{A21} \]

Step 20: Identify the Distributions

Step 20.1: Identify the Inverse-Gamma

The first factor in (A21) is:

\[ (\sigma^2)^{-\left(a_0 + \frac{n}{2} + 1\right)} \exp\left[-\frac{1}{\sigma^2}\left(b_0 + \frac{f(\theta, \bar{y})}{2}\right)\right] \]

This is proportional to an Inverse-Gamma density with:

• Shape: \(a_1 = a_0 + \frac{n}{2}\)
• Scale: \(b_1 = b_0 + \frac{f(\theta, \bar{y})}{2}\)

Therefore:

\[ \text{First factor} \propto IG(\sigma^2 \mid a_1, b_1) \tag{A22} \]

Step 20.2: Identify the Normal

The second factor in (A21) is:

\[ (\sigma^2)^{-1/2} \exp\left[-\frac{n_0 + n}{2\sigma^2}\left(\mu - \frac{n_0\theta + n\bar{y}}{n_0 + n}\right)^2\right] \]

Recall the Normal density:

\[ N\left(\mu \,\middle|\, m, \frac{\sigma^2}{k}\right) \propto (\sigma^2)^{-1/2} \exp\left[-\frac{k}{2\sigma^2}(\mu - m)^2\right] \]

Here, \(k = n_0 + n\) and \(m = \frac{n_0\theta + n\bar{y}}{n_0 + n}\).

Therefore:

\[ \text{Second factor} \propto N\left(\mu \,\middle|\, \frac{n_0\theta + n\bar{y}}{n_0 + n}, \frac{\sigma^2}{n_0 + n}\right) \tag{A23} \]

---

Part 6: The Final Factorization

Step 21: Final Result

Combining (A21), (A22), and (A23):

\[ \boxed{ p(\mu, \sigma^2 \mid y) \propto \underbrace{IG(\sigma^2 \mid a_1, b_1)}_{p(\sigma^2 \mid y)} \times \underbrace{N\left(\mu \,\middle|\, \frac{n_0\theta + n\bar{y}}{n_0 + n}, \frac{\sigma^2}{n_0 + n}\right)}_{p(\mu \mid \sigma^2, y)} } \tag{A24} \]

where:

\[ \boxed{ a_1 = a_0 + \frac{n}{2}, \quad b_1 = b_0 + \frac{f(\theta, \bar{y})}{2} } \tag{A25} \]

and:

\[ \boxed{ f(\theta, \bar{y}) = (n-1)s^2 + \frac{n_0 n}{n_0 + n}(\bar{y} - \theta)^2 } \tag{A26} \]

---

Part 7: Interpretation of the Factorization

What This Factorization Reveals

1. The conditional posterior \(p(\mu \mid \sigma^2, y)\) is Normal with:

   • Mean: \(\frac{n_0\theta + n\bar{y}}{n_0 + n}\) (a weighted average of prior mean \(\theta\) and sample mean \(\bar{y}\))
   • Variance: \(\frac{\sigma^2}{n_0 + n}\) (still depends on \(\sigma^2\))

2. The marginal posterior \(p(\sigma^2 \mid y)\) is Inverse-Gamma with:

   • Shape: \(a_1 = a_0 + \frac{n}{2}\)
   • Scale: \(b_1 = b_0 + \frac{f(\theta, \bar{y})}{2}\)

3. Crucially: \(\mu\) and \(\sigma^2\) remain dependent in the posterior because the variance of \(\mu \mid \sigma^2, y\) still contains \(\sigma^2\).

4. This factorization enables:

   • Gibbs sampling: Alternate between sampling \(\mu \mid \sigma^2, y\) and \(\sigma^2 \mid \mu, y\)
   • Closed-form marginal posteriors: Easy to compute posterior summaries

---

Summary of Key Steps

Step	Operation	Result
1	Apply Bayes’ theorem	\(p(\mu,\sigma^2 \mid y) \propto p(\mu,\sigma^2) p(y \mid \mu,\sigma^2)\)
2	Factor joint prior	\(= p(\mu \mid \sigma^2) p(\sigma^2) p(y \mid \mu,\sigma^2)\)
3-6	Write densities and drop constants	\(= (\sigma^2)^{-(a_0+1)} e^{-b_0/\sigma^2} \times (\sigma^2)^{-1/2} e^{-n_0(\mu-\theta)^2/(2\sigma^2)} \times (\sigma^2)^{-n/2} e^{-\sum(y_i-\mu)^2/(2\sigma^2)}\)
7-10	Combine powers and exponentials	\(= (\sigma^2)^{-(a_0+1+(n+1)/2)} \exp[-\frac{1}{\sigma^2}(b_0 + \frac{1}{2}Q)]\)
11-12	Define \(Q(\mu,\theta,y)\)	\(Q = n_0(\mu-\theta)^2 + \sum(y_i-\mu)^2\)
13-15	Expand and complete square	\(Q = (n_0+n)(\mu - m)^2 + f(\theta,\bar{y})\)
16-19	Substitute and separate	\(= [(\sigma^2)^{-(a_0+n/2+1)} e^{-(b_0+f/2)/\sigma^2}] \times [(\sigma^2)^{-1/2} e^{-(n_0+n)(\mu-m)^2/(2\sigma^2)}]\)
20-21	Identify distributions	\(= IG(\sigma^2 \mid a_1,b_1) \times N(\mu \mid m, \sigma^2/(n_0+n))\)

---

Key Takeaway

The factorization works because of the “completing the square” technique, which separates the \(\mu\) terms from the \(\sigma^2\) terms in the exponential. This reveals that the joint posterior is proportional to an Inverse-Gamma (for \(\sigma^2\)) times a Normal (for \(\mu\) given \(\sigma^2\)). This is a direct consequence of the conditional conjugacy built into the hierarchical prior.

---

–

Appendix: R Code for Gibbs Sampling

# Gibbs sampler for the hierarchical model
gibbs_sampler <- function(y, theta, n0, a0, b0, n_iter = 10000, burnin = 1000) {
  n <- length(y)
  ybar <- mean(y)
  
  # Initialize
  sigma2 <- 1
  mu <- ybar
  
  # Storage
  mu_samples <- numeric(n_iter)
  sigma2_samples <- numeric(n_iter)
  
  for (i in 1:n_iter) {
    # Sample mu | sigma2, y
    mu_mean <- (n0 * theta + n * ybar) / (n0 + n)
    mu_var <- sigma2 / (n0 + n)
    mu <- rnorm(1, mu_mean, sqrt(mu_var))
    
    # Sample sigma2 | mu, y
    a1 <- a0 + n / 2
    b1 <- b0 + 0.5 * (n0 * (mu - theta)^2 + sum((y - mu)^2))
    sigma2 <- 1 / rgamma(1, a1, b1)
    
    # Store
    mu_samples[i] <- mu
    sigma2_samples[i] <- sigma2
  }
  
  # Remove burnin
  mu_samples <- mu_samples[-(1:burnin)]
  sigma2_samples <- sigma2_samples[-(1:burnin)]
  
  return(list(mu = mu_samples, sigma2 = sigma2_samples))
}

# Example usage
set.seed(123)
n <- 30
theta <- 5
n0 <- 2
a0 <- 2
b0 <- 1
true_mu <- 3
true_sigma2 <- 4
y <- rnorm(n, true_mu, sqrt(true_sigma2))

samples <- gibbs_sampler(y, theta, n0, a0, b0, 10000)

# Posterior summaries
cat("Posterior mean of mu:", mean(samples$mu), "\n")

## Posterior mean of mu: 3.032456

cat("Posterior mean of sigma^2:", mean(samples$sigma2), "\n")

## Posterior mean of sigma^2: 3.919923

# Credible intervals
cat("95% CI for mu:", quantile(samples$mu, c(0.025, 0.975)), "\n")

## 95% CI for mu: 2.358915 3.719508

cat("95% CI for sigma^2:", quantile(samples$sigma2, c(0.025, 0.975)), "\n")

## 95% CI for sigma^2: 2.401608 6.376596

# Trace plots
par(mfrow = c(2, 2))
plot(samples$mu, type = "l", main = "Trace of mu", xlab = "Iteration", ylab = "mu")
hist(samples$mu, main = "Posterior of mu", breaks = 30, probability = TRUE)
plot(samples$sigma2, type = "l", main = "Trace of sigma^2", xlab = "Iteration", ylab = "sigma^2")
hist(samples$sigma2, main = "Posterior of sigma^2", breaks = 30, probability = TRUE)