Wishart Distribution and Bartlett’s Decomposition

The Wishart distribution is the multivariate generalization of chi-squared. It describes the distribution of sample covariance matrices. Bartlett’s decomposition provides an efficient method using simpler distributions.

Multivariate Normal Distribution

Definition

A random vector \(\mathbf{X} \in \mathbb{R}^p\) follows a multivariate normal distribution, written \(\mathbf{X} \sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), if its PDF is:

\[ f(\mathbf{x}) = \frac{1}{(2\pi)^{p/2} |\boldsymbol{\Sigma}|^{1/2}} \exp\left( -\frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right) \]

Parameters

Parameter	Dimension	Meaning
\(\boldsymbol{\mu}\)	\(p \times 1\)	Mean vector
\(\boldsymbol{\Sigma}\)	\(p \times p\)	Covariance matrix (must be PD [positive definite])

Key Properties

1. Linear transformation: If \(\mathbf{X} \sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), then:

\[ \mathbf{A}\mathbf{X} + \mathbf{b} \sim N_q(\mathbf{A}\boldsymbol{\mu} + \mathbf{b}, \mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^\top) \]

2. Marginals are normal: Any subset of variables is also multivariate normal

3. Conditionals are normal: \(\mathbf{X}_1 | \mathbf{X}_2 = \mathbf{x}_2\) is also normal

4. Mahalanobis distance: The quadratic form follows chi-squared:

\[ (\mathbf{X} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{X} - \boldsymbol{\mu}) \sim \chi^2_p \]

---

Wishart Distribution

Motivation

The Wishart distribution is the multivariate generalization of chi-squared. It describes the distribution of sample covariance matrices.

Univariate	Multivariate
\(\chi^2\) distribution	Wishart distribution
Sum of squared normals	Sum of outer products
Scalar	Matrix

Definition

Let \(\mathbf{X}_1, \ldots, \mathbf{X}_n \stackrel{iid}{\sim} N_p(\mathbf{0}, \boldsymbol{\Sigma})\). Then:

\[ \mathbf{W} = \sum_{i=1}^{n} \mathbf{X}_i \mathbf{X}_i^\top = \mathbf{X}^\top \mathbf{X} \sim W_p(\boldsymbol{\Sigma}, n) \]

where \(\mathbf{X}\) is the \(n \times p\) data matrix.

Parameters

Parameter	Meaning
\(\boldsymbol{\Sigma}\)	Scale matrix (\(p \times p\), PD)
\(n\)	Degrees of freedom (\(n \geq p\))

PDF

For \(\mathbf{W} \sim W_p(\boldsymbol{\Sigma}, n)\) with \(n \geq p\):

\[ f(\mathbf{W}) = \frac{|\mathbf{W}|^{(n-p-1)/2} \exp\left(-\frac{1}{2} \text{tr}(\boldsymbol{\Sigma}^{-1}\mathbf{W})\right)}{2^{np/2} |\boldsymbol{\Sigma}|^{n/2} \Gamma_p(n/2)} \]

where \(\Gamma_p\) is the multivariate gamma function.

Key Properties

1. Mean:

\[ E[\mathbf{W}] = n \boldsymbol{\Sigma} \]

2. Sample covariance connection: If \(\mathbf{S} = \frac{1}{n-1}\sum_{i=1}^{n}(\mathbf{X}_i - \bar{\mathbf{X}})(\mathbf{X}_i - \bar{\mathbf{X}})^\top\), then:

\[ (n-1)\mathbf{S} \sim W_p(\boldsymbol{\Sigma}, n-1) \]

3. Additivity: If \(\mathbf{W}_1 \sim W_p(\boldsymbol{\Sigma}, n_1)\) and \(\mathbf{W}_2 \sim W_p(\boldsymbol{\Sigma}, n_2)\) are independent:

\[ \mathbf{W}_1 + \mathbf{W}_2 \sim W_p(\boldsymbol{\Sigma}, n_1 + n_2) \]

4. Univariate case: \(W_1(\sigma^2, n) = \sigma^2 \chi^2_n\)

---

Connection Summary

\[ \boxed{ \begin{aligned} \text{Normal} &\xrightarrow{\text{square}} \chi^2 \\ \text{Multivariate Normal} &\xrightarrow{\text{outer product}} \text{Wishart} \end{aligned} } \]

Concept	Univariate	Multivariate
Distribution of data	\(N(\mu, \sigma^2)\)	\(N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)
Distribution of variance estimate	\(\chi^2\)	Wishart
Parameter	\(\sigma^2\) (scalar)	\(\boldsymbol{\Sigma}\) (matrix)

Multivariate Normal and Wishart in Bayesian Statistics

---

1. The Bayesian Setup

In Bayesian inference for multivariate normal data:

\[ \mathbf{X}_1, \ldots, \mathbf{X}_n \mid \boldsymbol{\mu}, \boldsymbol{\Sigma} \stackrel{iid}{\sim} N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \]

We need priors for:

• \(\boldsymbol{\mu}\) (mean vector)
• \(\boldsymbol{\Sigma}\) (covariance matrix) or \(\boldsymbol{\Omega} = \boldsymbol{\Sigma}^{-1}\) (precision matrix)

---

2. Conjugate Priors

Parameter	Conjugate Prior	Result
\(\boldsymbol{\mu}\) (known \(\boldsymbol{\Sigma}\))	Multivariate Normal	Normal posterior
\(\boldsymbol{\Sigma}\) (known \(\boldsymbol{\mu}\))	Inverse Wishart	Inverse Wishart posterior
\(\boldsymbol{\Omega}\) (known \(\boldsymbol{\mu}\))	Wishart	Wishart posterior
\((\boldsymbol{\mu}, \boldsymbol{\Sigma})\) jointly	Normal-Inverse-Wishart	Normal-Inverse-Wishart posterior

---

3. Case 1: Unknown μ, Known Σ

Prior:

\[ \boldsymbol{\mu} \sim N_p(\boldsymbol{\mu}_0, \boldsymbol{\Lambda}_0) \]

Posterior:

\[ \boldsymbol{\mu} \mid \mathbf{X}_1, \ldots, \mathbf{X}_n \sim N_p(\boldsymbol{\mu}_n, \boldsymbol{\Lambda}_n) \]

where:

\[ \boldsymbol{\Lambda}_n^{-1} = \boldsymbol{\Lambda}_0^{-1} + n\boldsymbol{\Sigma}^{-1} \]

\[ \boldsymbol{\mu}_n = \boldsymbol{\Lambda}_n \left( \boldsymbol{\Lambda}_0^{-1} \boldsymbol{\mu}_0 + n\boldsymbol{\Sigma}^{-1} \bar{\mathbf{X}} \right) \]

Interpretation: Posterior mean is a precision-weighted average of prior mean and sample mean.

---

4. Case 2: Unknown Σ, Known μ

4a. Prior on Covariance Matrix

Prior:

\[ \boldsymbol{\Sigma} \sim W_p^{-1}(\boldsymbol{\Psi}, \nu) \]

This is the Inverse Wishart with: - \(\boldsymbol{\Psi}\): scale matrix (\(p \times p\), PD) - \(\nu\): degrees of freedom (\(\nu > p - 1\))

PDF:

\[ f(\boldsymbol{\Sigma}) \propto |\boldsymbol{\Sigma}|^{-(\nu + p + 1)/2} \exp\left( -\frac{1}{2} \text{tr}(\boldsymbol{\Psi} \boldsymbol{\Sigma}^{-1}) \right) \]

Prior Mean:

\[ E[\boldsymbol{\Sigma}] = \frac{\boldsymbol{\Psi}}{\nu - p - 1} \quad (\text{for } \nu > p + 1) \]

Posterior:

\[ \boldsymbol{\Sigma} \mid \mathbf{X}_1, \ldots, \mathbf{X}_n \sim W_p^{-1}(\boldsymbol{\Psi}_n, \nu_n) \]

where:

\[ \nu_n = \nu + n \]

\[ \boldsymbol{\Psi}_n = \boldsymbol{\Psi} + \sum_{i=1}^{n} (\mathbf{X}_i - \boldsymbol{\mu})(\mathbf{X}_i - \boldsymbol{\mu})^\top \]

4b. Prior on Precision Matrix

Alternatively, place a Wishart prior directly on \(\boldsymbol{\Omega} = \boldsymbol{\Sigma}^{-1}\):

\[ \boldsymbol{\Omega} \sim W_p(\mathbf{V}, n_0) \]

Posterior:

\[ \boldsymbol{\Omega} \mid \mathbf{X}_1, \ldots, \mathbf{X}_n \sim W_p(\mathbf{V}_n, n_0 + n) \]

where:

\[ \mathbf{V}_n^{-1} = \mathbf{V}^{-1} + \sum_{i=1}^{n} (\mathbf{X}_i - \boldsymbol{\mu})(\mathbf{X}_i - \boldsymbol{\mu})^\top \]

---

5. Case 3: Both μ and Σ Unknown

Use the Normal-Inverse-Wishart (NIW) prior:

\[ (\boldsymbol{\mu}, \boldsymbol{\Sigma}) \sim \text{NIW}(\boldsymbol{\mu}_0, \kappa_0, \boldsymbol{\Psi}, \nu) \]

This is defined as:

\[ \boldsymbol{\Sigma} \sim W_p^{-1}(\boldsymbol{\Psi}, \nu) \]

\[ \boldsymbol{\mu} \mid \boldsymbol{\Sigma} \sim N_p\left(\boldsymbol{\mu}_0, \frac{\boldsymbol{\Sigma}}{\kappa_0}\right) \]

Hyperparameters:

Parameter	Interpretation
\(\boldsymbol{\mu}_0\)	Prior mean of \(\boldsymbol{\mu}\)
\(\kappa_0\)	“Prior sample size” for mean (strength of belief)
\(\boldsymbol{\Psi}\)	Prior scale matrix for \(\boldsymbol{\Sigma}\)
\(\nu\)	“Prior sample size” for covariance (degrees of freedom)

Posterior:

\[ (\boldsymbol{\mu}, \boldsymbol{\Sigma}) \mid \text{data} \sim \text{NIW}(\boldsymbol{\mu}_n, \kappa_n, \boldsymbol{\Psi}_n, \nu_n) \]

where:

\[ \kappa_n = \kappa_0 + n \]

\[ \nu_n = \nu + n \]

\[ \boldsymbol{\mu}_n = \frac{\kappa_0 \boldsymbol{\mu}_0 + n \bar{\mathbf{X}}}{\kappa_0 + n} \]

\[ \boldsymbol{\Psi}_n = \boldsymbol{\Psi} + \mathbf{S} + \frac{\kappa_0 n}{\kappa_0 + n}(\bar{\mathbf{X}} - \boldsymbol{\mu}_0)(\bar{\mathbf{X}} - \boldsymbol{\mu}_0)^\top \]

where \(\mathbf{S} = \sum_{i=1}^{n}(\mathbf{X}_i - \bar{\mathbf{X}})(\mathbf{X}_i - \bar{\mathbf{X}})^\top\) is the scatter matrix.

---

6. Marginal and Conditional Posteriors

From the NIW posterior, we can derive:

Marginal for Σ:

\[ \boldsymbol{\Sigma} \mid \text{data} \sim W_p^{-1}(\boldsymbol{\Psi}_n, \nu_n) \]

Marginal for μ:

\[ \boldsymbol{\mu} \mid \text{data} \sim t_{\nu_n - p + 1}\left(\boldsymbol{\mu}_n, \frac{\boldsymbol{\Psi}_n}{\kappa_n(\nu_n - p + 1)}\right) \]

This is a multivariate t-distribution (heavier tails than normal due to uncertainty in \(\boldsymbol{\Sigma}\)).

---

7. Choosing Priors

Weakly Informative / Non-informative

Setting	Effect
\(\kappa_0 \to 0\)	Vague prior on \(\boldsymbol{\mu}\)
\(\nu \to p - 1\)	Vague prior on \(\boldsymbol{\Sigma}\) (limit of validity)
\(\boldsymbol{\Psi} = \mathbf{I}\)	No prior preference for covariance structure

Jeffreys Prior

\[ p(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \propto |\boldsymbol{\Sigma}|^{-(p+1)/2} \]

This is improper but yields proper posteriors when \(n > p\).

---

8. Why Inverse Wishart for Σ?

Reason	Explanation
Conjugacy	Posterior has same form as prior
Guarantees PD	Samples from Inverse Wishart are always positive definite
Interpretable	\(\nu\) acts like prior sample size

Limitations: - Couples variance and correlation (can’t specify separately) - May be too restrictive for high-dimensional problems

Alternatives for high dimensions: - LKJ prior (for correlation matrix) - Decomposition: \(\boldsymbol{\Sigma} = \text{diag}(\boldsymbol{\sigma}) \mathbf{R} \text{diag}(\boldsymbol{\sigma})\)

---

9. R Example: Bayesian Inference with NIW

library(MCMCpack)  # For riwish (inverse Wishart)
library(MASS)

# Simulate data
set.seed(42)
n <- 50
p <- 2
mu_true <- c(2, -1)
Sigma_true <- matrix(c(1, 0.6, 0.6, 1), p, p)
X <- mvrnorm(n, mu_true, Sigma_true)

# Sufficient statistics
X_bar <- colMeans(X)
S <- cov(X) * (n - 1)  # Scatter matrix

# Prior hyperparameters (weakly informative)
mu_0 <- rep(0, p)
kappa_0 <- 0.1
nu <- p + 2
Psi <- diag(p)

# Posterior hyperparameters
kappa_n <- kappa_0 + n
nu_n <- nu + n
mu_n <- (kappa_0 * mu_0 + n * X_bar) / kappa_n
Psi_n <- Psi + S + (kappa_0 * n / kappa_n) * tcrossprod(X_bar - mu_0)

# Sample from posterior
n_samples <- 5000
Sigma_samples <- array(dim = c(p, p, n_samples))
mu_samples <- matrix(nrow = n_samples, ncol = p)

for (i in 1:n_samples) {
  # Sample Sigma from Inverse Wishart
  Sigma_samples[,,i] <- riwish(nu_n, Psi_n)
  # Sample mu from conditional Normal
  mu_samples[i,] <- mvrnorm(1, mu_n, Sigma_samples[,,i] / kappa_n)
}

# Posterior summaries
cat("True μ:", mu_true, "\n")
cat("Posterior mean of μ:", colMeans(mu_samples), "\n")
cat("\nTrue Σ:\n"); print(Sigma_true)
cat("\nPosterior mean of Σ:\n"); print(apply(Sigma_samples, 1:2, mean))

Bartlett’s Decomposition of Wishart Distribution

The Problem

Sampling directly from a Wishart or Inverse Wishart distribution is complex. Bartlett’s decomposition provides an efficient method using simpler distributions.

---

The Decomposition

If \(\mathbf{W} \sim W_p(\mathbf{I}, n)\) (Wishart with identity scale matrix), then:

\[ \mathbf{W} = \mathbf{L} \mathbf{L}^\top \]

where \(\mathbf{L}\) is a lower triangular matrix with:

\[ \mathbf{L} = \begin{pmatrix} \sqrt{c_1} & 0 & 0 & \cdots & 0 \\ z_{21} & \sqrt{c_2} & 0 & \cdots & 0 \\ z_{31} & z_{32} & \sqrt{c_3} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ z_{p1} & z_{p2} & z_{p3} & \cdots & \sqrt{c_p} \end{pmatrix} \]

Elements are independent with:

Element	Distribution
\(c_i\) (diagonal)	\(c_i \sim \chi^2_{n-i+1}\)
\(z_{ij}\) (off-diagonal, \(i > j\))	\(z_{ij} \sim N(0, 1)\)

---

General Scale Matrix

For \(\mathbf{W} \sim W_p(\boldsymbol{\Sigma}, n)\), let \(\boldsymbol{\Sigma} = \mathbf{C}\mathbf{C}^\top\) (Cholesky decomposition). Then:

\[ \mathbf{W} = \mathbf{C} \mathbf{L} \mathbf{L}^\top \mathbf{C}^\top \]

where \(\mathbf{L}\) is constructed as above.

---

Inverse Wishart Sampling

For \(\boldsymbol{\Sigma} \sim W_p^{-1}(\boldsymbol{\Psi}, \nu)\):

1. Compute Cholesky: \(\boldsymbol{\Psi}^{-1} = \mathbf{C}\mathbf{C}^\top\)
2. Generate \(\mathbf{L}\) using Bartlett’s decomposition with df = \(\nu\)
3. Compute \(\mathbf{A} = \mathbf{C}\mathbf{L}\)
4. Return \(\boldsymbol{\Sigma} = (\mathbf{A}\mathbf{A}^\top)^{-1}\)

---

R Implementation

# Sample from Wishart using Bartlett's decomposition
rwishart_bartlett <- function(n, df, Sigma) {
  p <- nrow(Sigma)
  
  # Cholesky of Sigma
  C <- chol(Sigma)
  
  # Construct L (lower triangular)
  L <- matrix(0, p, p)
  
  # Diagonal: sqrt of chi-squared
  for (i in 1:p) {
    L[i, i] <- sqrt(rchisq(1, df - i + 1))
  }
  
  # Off-diagonal: standard normal
  for (i in 2:p) {
    for (j in 1:(i-1)) {
      L[i, j] <- rnorm(1)
    }
  }
  
  # W = C' * L * L' * C = (L' * C)' * (L' * C)
  A <- t(L) %*% C
  W <- t(A) %*% A
  
  return(W)
}

# Sample from Inverse Wishart
rinvwishart_bartlett <- function(n, df, Psi) {
  W <- rwishart_bartlett(1, df, solve(Psi))
  return(solve(W))
}

# Test
set.seed(42)
Sigma <- matrix(c(1, 0.5, 0.5, 1), 2, 2)
df <- 10

# Generate samples
samples <- replicate(5000, rwishart_bartlett(1, df, Sigma), simplify = FALSE)
mean_W <- Reduce(`+`, samples) / length(samples)

cat("Expected E[W] = df * Sigma:\n")
print(df * Sigma)
cat("\nSample mean:\n")
print(round(mean_W, 2))

Why Bartlett’s Decomposition Matters

Advantage	Explanation
Efficiency	Only need \(\chi^2\) and \(N(0,1)\) samples
Numerical stability	Works through Cholesky factors
Guarantees PD	Product \(\mathbf{LL}^\top\) is always PD
MCMC friendly	Easy to implement in Gibbs samplers

Connection to Posterior

In Bayesian inference with NIW prior, the posterior for \(\boldsymbol{\Sigma}\) is:

\[ \boldsymbol{\Sigma} \mid \text{data} \sim W_p^{-1}(\boldsymbol{\Psi}_n, \nu_n) \]

Bartlett’s decomposition allows efficient sampling from this posterior in each MCMC iteration.

A New Method for Sampling from the Wishart Distribution

A new method for sampling from the Wishart distribution is introduced by representing the Wishart distributed random matrix as a function of independent multivariate normal-gamma random vectors. Building on this, an efficient monotone data augmentation (MDA) algorithm is developed for Bayesian multivariate linear regression.

For longitudinal outcomes, this method is easier to implement and interpret than that based on Bartlett’s decomposition. The proposed algorithm is both computationally efficient and practical, a benefit that is illustrated by the analysis of an antidepressant trial.

The decomposition methodology proposed by Tang (2015) was employed in the following R code to efficiently sample from the Wishart distribution by representing the random matrix as a function of independent multivariate normal and Gamma random vectors.

It must be noted that this implementation reflects my best effort to reproduce the method and may not be perfectly accurate.

###########################################
# USING R PACKAGE TO PRODUCE A "BENCHMARK"
###########################################
library(rWishart)

set.seed(12345)

rep=1e6

# Wishart random matrix using R package
A <- matrix(c(1,1,1,2), 2, 2, byrow = TRUE)
samplewishart1 <- rWishart(rep, 2, A)
approx1 <- apply(samplewishart1, 1:2, mean) # Check properties
cat("E[W] using R package: \n")
print(approx1)

###########################################
# Tang (2015)
###########################################
# Wishart random matrix using algorithm from Tang 2015
n0=2
p=2

# consider k=2
k=2

A <- matrix(c(2,-1,-1,1), 2, 2, byrow = TRUE)
A_inv <- solve(A)
L <- chol(A_inv)
L_t <- t(L)
L_inv <- solve(L)
print(L_t%*%L)
print(L_inv%*%t(L_inv))

a11 <- a12 <- a21 <- a22 <- 0

for (i in seq_len(rep)) {
  
  t11 <- L[1,1]*sqrt(rchisq(n = 1,df = n0-p+1))
  t22 <- sqrt(rchisq(n = 1,df = n0-p+2))
  t21 <- rnorm(1,mean=0, sd=1)
  
  T <- matrix(c(t11,0,t21,t22), 2, 2, byrow = TRUE)
  H <- T%*%L
  omega <- t(H)%*%H
  
  a11 <- a11 + omega[1, 1]
  a12 <- a12 + omega[1, 2]
  a21 <- a21 + omega[2, 1]
  a22 <- a22 + omega[2, 2]
}

avgM <- matrix(c(a11, a12, a21, a22) / rep, 2, 2, byrow = TRUE)
cat("E[W] using algorithm by Tang 2015: \n")
print(avgM) # Check properties