Problem Setup

Consider the conjugate Bayesian linear model:

Likelihood: \[y_i \mid \beta, \sigma^2 \stackrel{\text{ind}}{\sim} N(x_i^T \beta, \sigma^2), \quad i = 1, \ldots, n\]

In matrix form: \[y \mid \beta, \sigma^2 \sim N(X\beta, \sigma^2 I_n)\]

Prior: \[\beta \mid \sigma^2 \sim N(\mu, \sigma^2 M_0)\] where \(\mu = M_0 m_0\)

\[\sigma^2 \sim IG(a_0, b_0)\]

Goal: Derive the joint posterior:

\[p(\beta, \sigma^2 \mid y) = IG(\sigma^2 \mid a^*, b^*) \times N(\beta \mid Mm, \sigma^2 M)\]

where: - \(m = M_0^{-1}\mu + X^T y\) - \(M^{-1} = M_0^{-1} + X^T X\) - \(a^* = a_0 + \frac{n}{2}\) - \(b^* = b_0 + \frac{1}{2}(\mu^T M_0^{-1}\mu + y^T y - m^T M m)\)

Step-by-Step Derivation

Step 1: Write the Likelihood Function

The likelihood for all observations is:

\[p(y \mid \beta, \sigma^2) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(y_i - x_i^T \beta)^2}{2\sigma^2}\right\}\]

\[p(y \mid \beta, \sigma^2) = (2\pi\sigma^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - x_i^T \beta)^2\right\}\]

In matrix notation:

\[p(y \mid \beta, \sigma^2) = (2\pi\sigma^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma^2}(y - X\beta)^T(y - X\beta)\right\}\]

Step 2: Write the Prior Densities

Prior for \(\beta \mid \sigma^2\):

\[p(\beta \mid \sigma^2) = (2\pi\sigma^2)^{-p/2} |M_0|^{-1/2} \exp\left\{-\frac{1}{2\sigma^2}(\beta - \mu)^T M_0^{-1}(\beta - \mu)\right\}\]

where \(p\) is the dimension of \(\beta\)，it is not the pdf.

Prior for \(\sigma^2\):

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

Step 3: Apply Bayes’ Theorem

The joint posterior is proportional to likelihood × prior:

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

Substituting the expressions:

\[p(\beta, \sigma^2 \mid y) \propto (2\pi\sigma^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma^2}(y - X\beta)^T(y - X\beta)\right\}\]

\[\times (2\pi\sigma^2)^{-p/2} |M_0|^{-1/2} \exp\left\{-\frac{1}{2\sigma^2}(\beta - \mu)^T M_0^{-1}(\beta - \mu)\right\}\]

\[\times \frac{b_0^{a_0}}{\Gamma(a_0)} (\sigma^2)^{-a_0-1} \exp\left\{-\frac{b_0}{\sigma^2}\right\}\]

Step 4: Combine Terms

Collect all constants and combine powers of \(\sigma^2\):

Constants (not involving \(\beta\) or \(\sigma^2\)):

\[(2\pi)^{-(n+p)/2} |M_0|^{-1/2} \frac{b_0^{a_0}}{\Gamma(a_0)}\]

Power of \(\sigma^2\):

\[(\sigma^2)^{-n/2-p/2-a_0-1} = (\sigma^2)^{-(a_0 + n/2) - p/2 - 1}\]

Exponential term:

\[\exp\left\{-\frac{1}{\sigma^2}\left[b_0 + \frac{1}{2}(y - X\beta)^T(y - X\beta) + \frac{1}{2}(\beta - \mu)^T M_0^{-1}(\beta - \mu)\right]\right\}\]

Thus:

\[p(\beta, \sigma^2 \mid y) \propto (\sigma^2)^{-a^* - p/2 - 1} \exp\left\{-\frac{1}{\sigma^2}\left[b_0 + \frac{1}{2}Q(\beta)\right]\right\}\]

where \(a^* = a_0 + n/2\) and:

\[Q(\beta) = (y - X\beta)^T(y - X\beta) + (\beta - \mu)^T M_0^{-1}(\beta - \mu)\]

Step 5: Complete the Square in \(\beta\)

This is the key algebraic step. Expand \(Q(\beta)\):

\[Q(\beta) = (y - X\beta)^T(y - X\beta) + (\beta - \mu)^T M_0^{-1}(\beta - \mu)\]

Expand the first term:

\[(y - X\beta)^T(y - X\beta) = y^T y - 2\beta^T X^T y + \beta^T X^T X \beta\]

Expand the second term:

\[(\beta - \mu)^T M_0^{-1}(\beta - \mu) = \beta^T M_0^{-1}\beta - 2\beta^T M_0^{-1}\mu + \mu^T M_0^{-1}\mu\]

Combine:

\[Q(\beta) = y^T y - 2\beta^T X^T y + \beta^T X^T X\beta + \beta^T M_0^{-1}\beta - 2\beta^T M_0^{-1}\mu + \mu^T M_0^{-1}\mu\]

\[Q(\beta) = \beta^T(X^T X + M_0^{-1})\beta - 2\beta^T(X^T y + M_0^{-1}\mu) + y^T y + \mu^T M_0^{-1}\mu\]

Step 6: Define the Posterior Precision Matrix

Let:

\[M^{-1} = M_0^{-1} + X^T X\]

And define:

\[m = M_0^{-1}\mu + X^T y\]

Then:

\[Q(\beta) = \beta^T M^{-1}\beta - 2\beta^T m + y^T y + \mu^T M_0^{-1}\mu\]

Step 7: Complete the Square (Continued)

We want to express \(Q(\beta)\) in the form:

\[Q(\beta) = (\beta - Mm)^T M^{-1}(\beta - Mm) + \text{constant}\]

Expand the proposed form:

\[(\beta - Mm)^T M^{-1}(\beta - Mm) = \beta^T M^{-1}\beta - 2\beta^T M^{-1}Mm + m^T Mm\]

\[= \beta^T M^{-1}\beta - 2\beta^T m + m^T Mm\]

Compare with our expression:

\[Q(\beta) = \beta^T M^{-1}\beta - 2\beta^T m + y^T y + \mu^T M_0^{-1}\mu\]

Therefore, to match, we need:

\[y^T y + \mu^T M_0^{-1}\mu = m^T Mm + \text{constant}\]

So:

\[Q(\beta) = (\beta - Mm)^T M^{-1}(\beta - Mm) + y^T y + \mu^T M_0^{-1}\mu - m^T Mm\]

Step 8: Substitute Back into the Posterior

\[p(\beta, \sigma^2 \mid y) \propto (\sigma^2)^{-a^* - p/2 - 1}\]

\[\times \exp\left\{-\frac{1}{2\sigma^2}(\beta - Mm)^T M^{-1}(\beta - Mm)\right\}\]

\[\times \exp\left\{-\frac{1}{\sigma^2}\left[b_0 + \frac{1}{2}(y^T y + \mu^T M_0^{-1}\mu - m^T Mm)\right]\right\}\]

Step 9: Define \(b^*\)

Let:

\[b^* = b_0 + \frac{1}{2}(y^T y + \mu^T M_0^{-1}\mu - m^T Mm)\]

Then:

\[p(\beta, \sigma^2 \mid y) \propto (\sigma^2)^{-a^* - p/2 - 1}\]

\[\times \exp\left\{-\frac{1}{2\sigma^2}(\beta - Mm)^T M^{-1}(\beta - Mm)\right\}\]

\[\times \exp\left\{-\frac{b^*}{\sigma^2}\right\}\]

Step 10: Separate the Terms

For \(\sigma^2\):

\[p(\sigma^2 \mid y) \propto (\sigma^2)^{-a^* - 1} \exp\left\{-\frac{b^*}{\sigma^2}\right\}\]

This is the kernel of an Inverse Gamma distribution:

\[p(\sigma^2 \mid y) = IG(\sigma^2 \mid a^*, b^*)\]

For \(\beta \mid \sigma^2, y\):

\[p(\beta \mid \sigma^2, y) \propto (\sigma^2)^{-p/2} \exp\left\{-\frac{1}{2\sigma^2}(\beta - Mm)^T M^{-1}(\beta - Mm)\right\}\]

This is the kernel of a Normal distribution:

\[p(\beta \mid \sigma^2, y) = N(\beta \mid Mm, \sigma^2 M)\]

Step 11: Combine into the Final Form

\[p(\beta, \sigma^2 \mid y) = IG(\sigma^2 \mid a^*, b^*) \times N(\beta \mid Mm, \sigma^2 M)\]

where:

\[m = M_0^{-1}\mu + X^T y\]

\[M^{-1} = M_0^{-1} + X^T X\]

\[a^* = a_0 + \frac{n}{2}\]

\[b^* = b_0 + \frac{1}{2}(\mu^T M_0^{-1}\mu + y^T y - m^T Mm)\]

Summary Table

Quantity Expression
Posterior Precision Matrix \(M^{-1} = M_0^{-1} + X^T X\)
Posterior Mean Term \(m = M_0^{-1}\mu + X^T y\)
Posterior Mean of \(\beta\) \(E[\beta \mid \sigma^2, y] = M m\)
Posterior Shape \(a^* = a_0 + \frac{n}{2}\)
Posterior Rate \(b^* = b_0 + \frac{1}{2}(\mu^T M_0^{-1}\mu + y^T y - m^T Mm)\)

Key Insights

  1. The prior is conjugate — the posterior has the same Normal-Inverse Gamma form as the prior.

  2. The posterior precision \(M^{-1}\) is the sum of prior precision \(M_0^{-1}\) and data precision \(X^T X\).

  3. The term \(b^*\) updates the Inverse Gamma rate parameter by incorporating the residual sum of squares and prior information.

  4. The factorization shows that conditional on \(\sigma^2\), the posterior for \(\beta\) is Normal, and the marginal posterior for \(\sigma^2\) is Inverse Gamma.

R Code Implementation Example

# Example: Simulate data and compute posterior parameters

# Set seed for reproducibility
set.seed(123)

# Generate data
n <- 100
p <- 3
X <- matrix(rnorm(n*p), nrow = n, ncol = p)
beta_true <- c(1, 2, -1)
sigma2_true <- 2
y <- X %*% beta_true + rnorm(n, 0, sqrt(sigma2_true))

# Prior parameters
mu <- c(0, 0, 0)
M0 <- diag(10, p)
a0 <- 2
b0 <- 1

# Compute posterior parameters
M_inv <- solve(M0) + t(X) %*% X
M <- solve(M_inv)
m <- solve(M0) %*% mu + t(X) %*% y
a_star <- a0 + n/2
b_star <- b0 + 0.5 * (t(mu) %*% solve(M0) %*% mu + t(y) %*% y - t(m) %*% M %*% m)

# Display results
#Posterior Parameters:

print(M)
##              [,1]          [,2]          [,3]
## [1,] 0.0121912900  0.0006676702  0.0013305284
## [2,] 0.0006676702  0.0106971364 -0.0001013472
## [3,] 0.0013305284 -0.0001013472  0.0111520730
print(m)
##            [,1]
## [1,]   76.39042
## [2,]  187.38535
## [3,] -104.62023
cat("\na* =", a_star)
## 
## a* = 52
cat("\nb* =", b_star)
## 
## b* = 107.4956

Sampling from the Posterior

# Load required library
library(MASS)

# Sample from the posterior
n_samples <- 1000

# First sample sigma^2 from Inverse Gamma
sigma2_samples <- 1 / rgamma(n_samples, shape = a_star, rate = b_star)

# Then sample beta conditional on sigma^2
beta_samples <- matrix(NA, nrow = n_samples, ncol = p)
for(i in 1:n_samples) {
  beta_samples[i, ] <- mvrnorm(1, mu = M %*% m, Sigma = sigma2_samples[i] * M)
}

# Summary of posterior samples
cat("\nPosterior Summary:\n")
## 
## Posterior Summary:
cat("Beta posterior means:\n")
## Beta posterior means:
print(colMeans(beta_samples))
## [1]  0.9146038  2.0760449 -1.0851320
cat("\nBeta posterior standard deviations:\n")
## 
## Beta posterior standard deviations:
print(apply(beta_samples, 2, sd))
## [1] 0.1559646 0.1493832 0.1531620
cat("\nSigma^2 posterior mean:", mean(sigma2_samples))
## 
## Sigma^2 posterior mean: 2.115473
cat("\nSigma^2 posterior sd:", sd(sigma2_samples))
## 
## Sigma^2 posterior sd: 0.2962869
# Plot posterior distributions
par(mfrow = c(2, 2))
for(j in 1:p) {
  hist(beta_samples[, j], main = paste("Beta", j), xlab = "", breaks = 30)
  abline(v = beta_true[j], col = "red", lwd = 2)
}
hist(sigma2_samples, main = "Sigma^2", xlab = "", breaks = 30)
abline(v = sigma2_true, col = "red", lwd = 2)