Gibbs Sampling for Posterior Inference

• Gibbs sampling is a Markov chain Monte Carlo method that leverages conditional independence to sample from high-dimensional posterior distributions in Bayesian inference.
• It employs block updates that alternate between multivariate normal and gamma mixtures to efficiently navigate the parameter space in hierarchical models.
• Its geometric ergodicity guarantees rapid convergence and reliable error estimation, underpinning practical applications in rigorous posterior analysis.

Gibbs sampling is a Markov chain Monte Carlo (MCMC) algorithm that generates dependent samples from a high-dimensional probability distribution, typically the posterior distribution arising in Bayesian inference. For posterior inference in hierarchical models, Gibbs sampling exploits conditional independence structure to iteratively sample from low-dimensional (often univariate or blockwise) conditional distributions, enabling efficient exploration of the potentially complex and unimodal or multimodal target distribution. Gibbs sampling has become a foundational tool in Bayesian computation, as it is applicable to a wide range of posterior structures where direct sampling is infeasible.

1. Hierarchical General Linear Model and Posterior Structure

For the Bayesian hierarchical general linear model, the data vector $Y \in \mathbb{R}^n$ is modeled as

$$Y \mid \beta, u, \lambda_R, \lambda_D \sim N_n(X\beta + Zu, \lambda_R^{-1} I_n)$$

where $\beta$ are fixed effects, $u$ are random effects, $X$ and $Z$ are known design matrices, and $\lambda_R, \lambda_D$ are precision parameters (inverse variances) for the residual and random effects, respectively. Priors are placed on $\beta$, $u$, and the precisions, often as finite mixtures of normal and gamma densities.

The joint posterior takes the form

$$\pi(\xi, \lambda \mid y) \propto f(y \mid \xi, \lambda) f(\xi \mid \lambda) f(\lambda)\,,$$

with $\xi = [u^\top, \beta^\top]^\top$ and $\lambda = [\lambda_R, \lambda_D]^\top$. Due to the hierarchical dependence, this joint cannot be sampled directly, motivating block Gibbs updates.

2. Block Gibbs Sampler: Algorithmic Construction

The core Gibbs sampler alternates between updating the “location” block (parameters $\xi$) and the “precision” block ($\lambda$), each conditional on current values of the other. The full conditional for the precisions is a mixture of two independent gamma densities (possibly weighted), depending on quadratic forms of the current $\xi$: $$v_1(\xi) = (y - X\beta - Zu)^\top (y - X\beta - Zu)\,, \quad v_2(\xi) = u^\top u$$ The conditional density for the precisions is

$$f(\lambda \mid \xi, y) = \sum_{j=1}^s \sum_{l=1}^t \phi_j \psi_l f_{1j}(\lambda_R \mid \xi, y) f_{2l}(\lambda_D \mid \xi, y)$$

where $f_{1j}, f_{2l}$ are gamma densities parametrized by $(v_1(\xi), v_2(\xi))$. The full conditional for $\xi$ is a mixture of multivariate normal densities: $$\xi \mid \lambda, y \sim \sum_{i=1}^r \eta_i N_{k+p}(m_i, \Sigma^{-1})$$ with block diagonal precision and location determined by the current precisions and the prior parameters.

The algorithm repeatedly cycles through:

• $\xi$‐step: Sample $\xi$ from its multivariate normal mixture conditional on $\lambda$.
• $\lambda$‐step: Sample $\lambda$ from the independent gamma mixture conditional on $\xi$.

Two natural update orders are considered ($\Phi_1$: $\xi$ then $\lambda$, $\Phi_2$: $\lambda$ then $\xi$). Both yield Markov chains with $\pi(\xi, \lambda \mid y)$ as invariant distribution.

3. Geometric Ergodicity: Convergence and Theoretical Guarantees

The paper establishes that under mild and explicit conditions, the block Gibbs sampler is geometrically ergodic: the total variation distance to stationarity decreases at an exponential rate,

$$\|\mathcal{L}(X^n) - \pi\|_{\mathrm{TV}} \leq C \gamma^n\,,\quad 0 < \gamma < 1$$

A Lyapunov (drift) function $V(\xi) = v_1(\xi) + v_2(\xi)$ is constructed, and a drift condition is verified: $$E[V(\xi_{n+1}) \mid \xi_n] \leq \gamma V(\xi_n) + L$$ with explicit constants in terms of the dimensions and Gamma prior parameters. Geometric ergodicity is essential for two reasons:

• It implies a Markov chain central limit theorem (CLT) for ergodic averages.
• It underpins strong laws for variance estimation (justifying standard error calculations) and the validity of sequential stopping rules.

The proof employs minorization conditions to ensure the Markov chain can “jump” into low-drift regions regardless of current state.

4. Central Limit Theorem and Consistency of Error Estimates

If the chain is geometrically ergodic and the function of interest $g$ has finite $(2+\varepsilon)$ moments under $\pi$, the central limit theorem holds: $$\sqrt{n}\,(\hat{g}_n - E_\pi g) \xrightarrow{d} N(0, \sigma_g^2)$$ where

$$\hat{g}_n = \frac{1}{n} \sum_{i=0}^{n-1} g(\xi_i, \lambda_i)$$

The asymptotic variance $\sigma_g^2$ is generally nontrivial because the Markov chain samples are autocorrelated. The batch means estimator divides the samples into $a_n$ blocks of size $b_n$, computes block averages, and uses their empirical variance as a consistent estimator: $$\hat{\sigma}_n^2 = \frac{b_n}{a_n - 1} \sum_{j=1}^{a_n} \left(\bar{Y}_j - \bar{g}_n\right)^2$$ This estimator is shown to be consistent under the established geometric ergodicity and moment conditions, enabling practitioners to construct asymptotically valid confidence intervals,

$$\bar{g}_n \pm t_{a_n-1} \frac{\hat{\sigma}_n}{\sqrt{n}}$$

5. Implementation in Practice: Health Plan Cost Example

The approach is instantiated on US government health maintenance organization (HMO) cost data, with a Bayesian linear regression: $$y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \varepsilon_i\,, \quad \varepsilon_i \sim N(0, \lambda_R^{-1})$$ with hierarchical prior

$$y \mid \beta, \lambda_R \sim N_N(X\beta, \lambda_R^{-1}I_N),\, \beta \mid \lambda_R \sim N_3(b,B^{-1}),\, \lambda_R \sim \mathrm{Gamma}(r_1,r_2)$$

Only $\beta$ and $\lambda_R$ are sampled, alternating normal and gamma full conditional draws. The empirical Bayes approach is used to select hyperparameters based on least squares estimates.

Convergence diagnostics are based on monitoring Monte Carlo standard errors (via the batch means estimator) for posterior means of interest. In the reported implementation, a total of 16,831 iterations bring all confidence interval half-widths below pre-set thresholds, and precise posterior means with their MCSEs are reported.

6. Practical Impact and Summary

The block Gibbs sampler constructed for the Bayesian hierarchical general linear model and its extensions is validated both theoretically and practically:

• Correct and efficient implementation: Alternating between mixtures of multivariate normal and independent gamma full conditionals.
• Geometric ergodicity: Demonstrated by explicit drift functions and minorization, guaranteeing rapid convergence and enabling CLT for averages.
• Variance estimation via batch means: Under geometric ergodicity, the batch means estimator is consistent, supporting practical error quantification.
• Empirical validation: On real data, geometric ergodicity and automatic error estimates enable rigorous posterior inference using Markov chain samples.
• General applicability: The methodology is broadly applicable to many hierarchical Bayesian models with normal likelihoods and general mixture priors on regression and variance parameters, and is not tied to conjugacy. Theoretical results are provided in explicit, verifiable terms with direct implications for applied practice.

This approach equips practitioners with a methodologically solid framework for confident inference in models too complex for direct sampling, ensuring that inferences grounded on Gibbs sampler output can be assessed with the same rigor as if based on independent posterior samples (0712.3056).