Gibbs Sampler

Definition

The Gibbs sampler is a Markov Chain Monte Carlo (MCMC) algorithm that generates samples from a high-dimensional joint posterior $p(\theta_1, \ldots, \theta_n \mid Y)$ by iteratively drawing each block $\theta_i$ from its full conditional $p(\theta_i \mid \theta_{-i}, Y)$ , cycling until the chain converges to the target distribution.

Key Ideas

The algorithm requires only the full conditional posteriors — it never evaluates the joint posterior directly.
By the ergodic theorem, time averages over the chain converge almost surely to population expectations: $k^{-1}\sum_j g(x^{(j)}) \xrightarrow{\text{a.s.}} E[g(X)]$ .
Convergence rate depends heavily on posterior correlations between blocks: high correlation $\Rightarrow$ high autocorrelation in the chain $\Rightarrow$ many draws needed.
Variants (Griddy-Gibbs, blocked Gibbs, Metropolis-within-Gibbs) handle cases where full conditionals are non-standard or highly correlated.

How It Works

The Gibbs Sampling Algorithm

Given a joint density $f(x_1, \ldots, x_n)$ with known full conditional densities $f(x_i \mid x_{-i})$ :

Initialize $x^{(0)} = (x_1^{(0)}, \ldots, x_n^{(0)})$ arbitrarily.
At iteration $j$ , draw sequentially: $x_1^{(j)} \sim f(x_1 \mid x_2^{(j-1)}, \ldots, x_n^{(j-1)})$ $x_2^{(j)} \sim f(x_2 \mid x_1^{(j)}, x_3^{(j-1)}, \ldots, x_n^{(j-1)})$ $\vdots$ $x_n^{(j)} \sim f(x_n \mid x_1^{(j)}, \ldots, x_{n-1}^{(j)})$
Repeat for $j = 1, \ldots, k$ . Discard a burn-in period; use $\{x^{(j)}\}$ as approximate draws from $f(x)$ .

The resulting sequence is a Markov chain whose equilibrium distribution is $f(x)$ under mild regularity conditions (Geman and Geman 1984; Tierney 1994).

Why the Gibbs Sampler Works: Markov Chain Foundation (Casella-George 1992)

That the Gibbs sequence converges to the target marginal is a consequence of Markov chain theory (Casella and George 1992). In the discrete 2-variable case, define:

$A_{xy}$ : matrix of $P(X = x \mid Y = y)$ — transitions from $y$ -states to $x$ -states
$A_{yx}$ : matrix of $P(Y = y \mid X = x)$ — transitions from $x$ -states to $y$ -states

The $X$ -only subsequence $X_0, X_1, \ldots$ is a Markov chain with transition matrix: $A_{xx} = A_{xy} \cdot A_{yx}$

Provided all entries of $A_{xx}$ are positive (irreducibility), the chain converges to a unique stationary distribution $\mathbf{f}$ satisfying $\mathbf{f}\,A_{xx} = \mathbf{f}$ . One verifies directly that this $\mathbf{f}$ is the marginal distribution of $X$ — the iteration has nowhere else to go. In the continuous case, the same argument extends via the integral analogue (their eq. 3.5).

Conditionals determine marginals. If the joint $f(x,y)$ exists, then $f_X(x)$ satisfies the fixed-point integral equation: $f_X(x) = \int h(x,t)\, f_X(t)\, dt, \qquad h(x,t) = \int f(x|y)\, f(y|t)\, dy \tag{FP}$ The Gibbs iteration is the computational device for solving (FP).

Failure mode 1: improper marginals. When $f(x|y)$ and $f(y|x)$ do not correspond to any proper joint density, the Gibbs sequence fails to converge. Canonical example: $f(x|y) = ye^{-xy}$ , $f(y|x) = xe^{-xy}$ on $(0,\infty)^2$ . Equation (FP) yields $f_X(x) = 1/x$ , which is not a density. The sufficient condition for convergence is that $f_X$ is proper ( $\int f_X\,dx < \infty$ ); restricting the conditionals to a compact interval guarantees this. Formal conditions are analyzed by Hobert and Casella (1996).

Failure mode 2: proper conditionals, improper joint (Besag, Scollnik). A subtler failure: a set of proper full conditional distributions need not define a proper joint. Example (Gelfand 1994, cited by Scollnik): $f(x_1, x_2) \propto \exp(-(x_1+x_2)^2/2)$ on $\mathbb{R}^2$ gives two valid univariate Normal conditionals — $f(x_1|x_2) = \mathcal{N}(-x_2,1)$ and $f(x_2|x_1) = \mathcal{N}(-x_1,1)$ — but the joint is improper. Running Gibbs on these conditionals will either fail to converge or converge to a state without a valid probabilistic interpretation. The Besag (1974) result guarantees uniqueness only when the joint distribution is known to be proper; it does not certify existence from conditionals alone.

Gibbs Samplers in VAR Models

The Gibbs sampler appears in several distinct vector autoregression (VAR) contexts, each exploiting different closed-form conditional structures:

Waggoner-Zha (2000, 2003b) — Structural VAR. Partitions structural parameters into equation-by-equation blocks $(b_i, g_i)$ : the lag-coefficient block $g_i$ has a Gaussian conditional; the contemporaneous block $b_i$ has a conditional that decomposes into a Univariate Wishart (UW) draw and Normal draws via an orthonormal basis (LU decomposition). See Structural Identification and the detailed derivation below.

Villani (2008) — Steady state VAR. Three-block sampler over $(\Sigma, \Pi, \Psi)$ where all full conditionals are Inverse-Wishart or Normal. The $\Psi$ block update stabilizes the sampler near unit roots. See Steady State VAR.

Bauwens-Giot (1997) — Cointegrated VAR. Full conditionals for the cointegrating vector elements have no closed form; the Griddy-Gibbs algorithm is used (see below). See Cointegration.

Albert-Chib (1993) — Markov-switching AR. Data augmentation: the latent state sequence $S_T$ is added as a Gibbs block, converting Hamilton's 2^T mixture likelihood into tractable conditionals. States drawn by a single backward binary pass; regression and variance blocks are truncated Normal and Inverse-Gamma; transition probabilities are Beta. Identified by $\gamma > 0$ , $\omega > 0$ . Foundational univariate Bayesian precursor to the MS-VAR literature. See Markov-Switching VAR.

Chib (1996) — Joint state sampling for Markov mixture models. Replaces the $n$ -block single-move sampler of Albert-Chib (1993) with a joint draw of the full state sequence $S_n = (s_1, \ldots, s_n)$ in one Gibbs block. Forward pass: compute $n \times m$ filter matrix $F$ via $F_t \propto (F_{t-1}P) \odot d_t$ ( $d_t$ = row of component densities at $t$ ; $\odot$ element-wise). Backward pass: draw $s_n$ from last row of $F$ ; for $t = n-1, \ldots, 1$ draw $s_t \propto p(s_t|Y_t,\theta) \cdot P_{s_t,s_{t+1}}$ (eq. 7). Transition matrix $P$ : Dirichlet conjugate from counts $n_{ij}$ (eq. 10). Discrete-state counterpart of the Carter-Kohn (1994) simulation smoother; became the standard for Bayesian Markov-switching estimation. See Markov-Switching VAR and Chib (1996).

Albert-Chib (1993b) — Binary and polychotomous probit. Data augmentation for binary regression: introduce latent $Z_i \sim \mathcal{N}(x_i'\beta, 1)$ with $Y_i = \mathbf{1}[Z_i > 0]$ . Two-block Gibbs: $\beta|Z \sim \mathcal{N}_k((X'X)^{-1}X'Z,\,(X'X)^{-1})$ (eq. 5) and $Z_i|y_i,\beta$ truncated Normal at 0 (eq. 6). Extensions: t-link via four-block scale-mixture sampler (eqs. 7–11), hierarchical shrinkage prior (eqs. 13–14), ordered multinomial with uniform threshold block (eqs. 16–18), unordered multinomial probit via accept-reject multivariate truncated Normal (eqs. 19–21). See Binary Probit.

Chib-Greenberg (1998) — Multivariate probit for correlated binary data. The $J$ -variate response $Y_i=(Y_{i1},\ldots,Y_{iJ})'$ with $Z_i \sim N_J(X_i\beta,\Sigma)$ and $Y_{ij}=I(z_{ij}>0)$ ; $\Sigma$ must be in correlation form (unit diagonal) for identifiability — $p=J(J-1)/2$ free off-diagonal elements $\sigma$ . Three-block Gibbs: (1) $Z_i|y_i,\beta,\Sigma$ — MVN truncated to $B_i$ (product of half-lines from $y_i$ ), sampled by cycling through $J$ univariate truncated normals (Geweke 1991); (2) $\beta|Z,\Sigma$ — conjugate Normal (eq. 6): $\hat\beta=B^{-1}(B_0\beta_0+\sum X_i'\Sigma^{-1}Z_i)$ ; (3) $\sigma|y,Z,\beta$ — Metropolis-Hastings (M-H) with tailored independence/reflection chain (eq. 7): proposal mode from Newton-Raphson on $\log g(\sigma|Z,\beta)$ , covariance from negative inverse Hessian; for high-dimensional $\sigma$ , partition into blocks applied sequentially. Likelihood evaluation via the GHK (Geweke-Hajivassiliou-Keane) method (Cholesky $\Sigma=LL'$ , change of variables, recursive truncated-normal draws, eq. 11). Posterior ordinate $\hat\pi(\sigma^*|y)$ via kernel density (eq. 10) — this is the precursor approach that Chib-Jeliazkov (2001) superseded with the two-expectation M-H identity. See Chib-Greenberg (1998), Binary Probit, and Multinomial Probit.

Kadiyala-Karlsson (1997) — Normal-Diffuse and Extended Natural Conjugate (ENC) Bayesian VAR (BVAR). Two new Gibbs samplers derived for the Normal-Diffuse (2-block) and ENC (equation-by-equation) posteriors; compared against importance sampling across 3- and 7-variable models. Gibbs achieves relative numerical efficiency (RNE) ≈ 1 at all model sizes; importance sampling collapses for the 7-variable model. See Noninformative Prior for VAR and Vector Autoregression.

Gibbs Simulator for Restricted Structural VARs (Waggoner-Zha 2000)

Write the structural VAR as $y_t' A = x_t' F + \varepsilon_t'$ . Linear zero restrictions take the form $Q_i a_i = 0$ on contemporaneous coefficients and $R_i f_i = 0$ on lag coefficients, where $a_i$ and $f_i$ are the $i$ th columns of $A$ and $F$ . The free parameters are $a_i = U_i b_i$ and $f_i = V_i g_i$ , where $U_i \in \mathbb{R}^{n \times k_i}$ and $V_i$ are orthonormal bases for the null spaces of $Q_i$ and $R_i$ .

Joint posterior (Theorem 1). Under the Sims-Zha (1998) reference prior, the joint posterior factorizes as $p(b_1,\ldots,b_n) \times \prod_i p(g_i | b_i)$ where:

$p(b_1,\ldots,b_n) \propto \left|\det[U_1 b_1 \mid \cdots \mid U_n b_n]\right|^T \exp\!\left(-\frac{T}{2}\sum_i b_i' S_i^{-1} b_i\right)$

The determinant factor renders the joint distribution non-Gaussian — it lies on a curved ridge in parameter space. Importance sampling with a Gaussian proposal concentrates almost all weight on a single draw (Figure 3b: single spike), making it useless for simultaneous-equations models.

Theorem 2 (Gibbs step for $b_i$ ). Conditional on all other equations, $b_i = \sum_{j=1}^{k_i} \beta_j w_j$ where $\{w_j\}$ is an orthonormal basis in $\mathbb{R}^{k_i}$ computed at each step via LU decomposition of the current draw of the other equations' free parameters. The scalar coefficients distribute as:

$\beta_1 \sim \text{UW}[T^{-1},\, k_i + 1], \qquad \beta_j \sim \mathcal{N}(0,\, T^{-1}) \quad (j \geq 2)$

Univariate Wishart (UW) distribution. The density is:

$f(x) \propto x^{(k-1)/2} \exp\!\left(-\frac{T}{2} x\right), \quad x > 0$

A draw is obtained as $x = z^2$ where $z \sim \mathcal{N}(0, T^{-1})$ — a one-line computation.

Lag-coefficient block. Conditional on the current $b_i$ , the free lag parameters $g_i$ have a Gaussian conditional — standard ordinary least squares (OLS) on the restricted system — and are drawn directly.

Theorem 3 / Corollary 1 (Independence conditions). The Gibbs draws across equations are exactly independent iff there exist a permutation matrix $\Gamma$ and orthogonal matrix $\Pi$ such that $\Gamma A \Pi$ is upper triangular for all admissible $A$ . Block-recursive exclusion restrictions are sufficient: if the contemporaneous matrix $A$ is block-lower-triangular after permutation, all cross-block draws are independent. In this case there is no need to run long chains — effective sample size equals the number of draws.

Convergence validation. Applied to Sims's (1986) 6-variable VAR ( $R$ , $M1$ , $Y$ , $P$ , $U$ , $ID$ ; 1948:1–1989:3; 20 free structural parameters in $A$ ). Five chains started from overdispersed initial draws; Gelman-Rubin $\hat{R} \to 1.00$ for all parameters within ~500 sweeps. The monetary-policy and interest-rate-defense equations satisfy Corollary 1 and yield exactly independent draws.

Normal-Diffuse and ENC Gibbs for BVAR (Kadiyala-Karlsson 1997)

Normal-Diffuse 2-Block Gibbs

The prior $\gamma \sim N(\tilde\gamma, \tilde\Sigma)$ , $p(\Psi) \propto |\Psi|^{-(q+1)/2}$ yields two closed-form full conditionals, enabling direct Gibbs sampling without importance functions:

Block 1 — $\gamma \mid \Psi, Y$ : Normal with precision $\tilde\Sigma^{-1} + \Psi^{-1} \otimes Z'Z$ : $\gamma|\Psi, Y \sim N\!\left(\bar\gamma,\; (\tilde\Sigma^{-1} + \Psi^{-1} \otimes Z'Z)^{-1}\right)$ $\bar\gamma = (\tilde\Sigma^{-1} + \Psi^{-1} \otimes Z'Z)^{-1}\!\left(\tilde\Sigma^{-1}\tilde\gamma + (\Psi^{-1} \otimes Z')y\right), \quad y = \text{vec}(Y)$

Block 2 — $\Psi^{-1} \mid \gamma, Y$ : Wishart with scale $S(\Phi)^{-1}$ and $T$ degrees of freedom: $\Psi^{-1}|\gamma, Y \sim W\!\left(S(\Phi)^{-1},\; T\right), \quad S(\Phi) = (Y - Z\Phi)'(Y - Z\Phi)$

Both blocks are standard random-variate draws. The dominant cost is the Cholesky decomposition of the $mk \times mk$ precision matrix in Block 1: $O((mk)^3)$ per Gibbs step.

ENC Equation-by-Equation Gibbs

Under the Drèze-Morales (1976) ENC reparametrization (see Noninformative Prior for VAR), there is no standard-form block decomposition for the full posterior. However, each equation's full conditional $\gamma_i|\gamma_{\neq i}, Y$ is a multivariate $t$ by Drèze-Richard (1983) Lemma 6.4: $\gamma_i | \gamma_{\neq i}, Y \sim t_k\!\left(\mathbf{d}_i,\; \frac{q_{ii}}{T+\alpha-k}\,\mathbf{P}_{ii}^{-1},\; T+\alpha-k\right)$ where $\mathbf{d}_i$ , $q_{ii}$ , $\mathbf{P}_{ii}$ are equation-specific quantities derived from the ENC reparametrization and the current draw of all other equations. One full cycle iterates over all $m$ equations at cost $O(m(mk)^2)$ — more efficient than the Normal-Diffuse 2-block sampler when $m$ and $k$ are both large.

Relative Numerical Efficiency and Importance Sampling Collapse

The relative numerical efficiency of method $M$ is: $\text{RNE}(M) = \frac{\text{Var}_p[g(\theta)]}{\sigma^2_{g,M}}$ It equals the fraction of direct i.i.d. draws that would achieve the same estimator variance. Gibbs sampling achieves RNE ≈ 1 at all model sizes. Importance sampling achieves moderate RNE for small models (3 variables) but collapses to RNE ≈ 0 for the 7-variable model because condition (15) — requiring $T - mk \geq 0$ so that the IS normalizing constant is numerically stable — fails severely ( $T - mk \approx -159$ in the reported experiment).

Antithetic variates (GAIN): Drawing $z^{*(i)} = -z^{(i)}$ (sign-flipped Normal variate through the posterior mode) reduces IS variance: $\text{GAIN} = \frac{\text{Var}(\bar g_{2n})}{\text{Var}(\bar g_{a,n})} = \frac{1}{1+\rho}, \quad \rho = \text{Corr}\!\left(g(\theta^{(i)}),\, g(\theta^{*(i)})\right)$ For near-symmetric posteriors $\rho \approx -1$ and $\text{GAIN} \gg 1$ . Draws are obtained from the antithetic 2-0 poly- $t$ importance function via the Bauwens-Richard (1982) decomposition. See also Geweke (1988) for the general principle.

CPU scaling across methods:

Method	Cost per draw/step	7-var / 3-var overhead
Normal-Wishart / Diffuse	$O(mk)$	${\approx}4\times$
STFC importance function	$O((mk)^{3/2})$	${\approx}8\times$
ENC equation-by-equation Gibbs	$O(m(mk)^2)$	${\approx}30\times$
Normal-Diffuse 2-block Gibbs	$O((mk)^3)$	${\approx}140\times$

Full comparative experiment required >200 CPU hours (1997 hardware) per prior across 3- and 7-variable Canadian macro VARs.

Model-Augmented Gibbs Sampler (Carlin-Chib 1995)

When comparing $K$ models via MCMC, the natural approach of adding a discrete model indicator $M$ to the Gibbs scheme fails due to an absorbing state problem: at iteration $g$ , if $M^{(g)} = j$ forces parameter $\theta_i$ ( $i \neq j$ ) out of the active model, there is no mechanism to regenerate $\theta_i$ at the next step — the chain gets trapped, violating irreducibility.

Solution: pseudo-priors. Carlin and Chib (1995) augment the joint model with linking densities $p(\theta_j | M \neq j)$ for inactive parameters. Since $y \perp \theta_j | (M \neq j)$ , these densities are computationally transparent — they don't affect $p(y|M=j)$ — but they keep every parameter vector accessible at every step. The full conditional for the model indicator is:

$P(M=j \mid \theta, y) \propto f(y|\theta_j,M=j)\; p(\theta_j|M=j)\; \pi_j \tag{CC}$

This is a proper discrete distribution, computable at each iteration. All other blocks alternate between the true model posterior (when active) and the pseudo-prior (when inactive). Standard regularity conditions (Roberts and Smith 1993; Tierney 1994) then guarantee convergence.

Practical implementation.

Run $K$ individual model-specific Gibbs chains to approximate the model-specific posteriors $p(\theta_j|M=j, y)$ .
Fit first-order normal (or parametric) approximations to these posteriors; use them as pseudo-priors.
Run the combined sampler; estimate $\hat{P}(M=j|y) = G_j/G$ from visit counts.

Pseudo-prior accuracy determines mixing speed: well-matched pseudo-priors produce frequent moves between models; diffuse pseudo-priors lead to high autocorrelation in the $M$ -chain. Prior weights $\pi_j$ can be adjusted to correct imbalance (one model visited too rarely), without invalidating the Bayes factor estimate.

Relation to Marginal Data Density. The Carlin-Chib sampler targets the posterior model probabilities $P(M=j|y)$ directly; Chib (1995) instead estimates $\log p(y|M_j)$ from a single-model run via the identity $\log p(y) = \log L(\phi^*) + \log p(\phi^*) - \log p(\phi^*|y)$ . Both approaches give the same Bayes factors; Carlin-Chib is more direct for small $K$ , Chib scales better for large $K$ .

Carter-Kohn Multimove Algorithm (State-Space Models)

In VAR-TVP models, the latent state path $\beta_1, \ldots, \beta_T$ must be drawn jointly rather than one element at a time (single-move Gibbs has severe mixing problems). The multimove algorithm (Carter and Kohn 1994) proceeds in two passes:

Step 1 — Forward Kalman filter ( $t = 1, \ldots, T$ ): run the standard Kalman recursions to compute filtered means $\beta_{t|t}$ and covariances $Q_{t|t}$ conditional on $y_1, \ldots, y_t$ and the current draw of static parameters.

Step 2 — Backward simulation smoother ( $t = T, T-1, \ldots, 1$ ): draw sequentially from Gaussians whose moments are updated using the next-period state draw: $\beta_T \sim \mathcal{N}(\beta_{T|T},\, Q_{T|T})$ $\beta_t \mid \beta_{t+1} \sim \mathcal{N}\!\left(\beta_{t|t} + \phi Q_{t|t} Q_{t+1|t}^{-1}(\beta_{t+1} - \phi\beta_{t|t} - (1-\phi)\beta_0),\; Q_{t|t} - \phi^2 Q_{t|t} Q_{t+1|t}^{-1} Q_{t|t}\right)$

The draws from this pass are exact samples from the joint posterior $p(\beta_1, \ldots, \beta_T \mid y_T, \theta)$ . The bottleneck is the $K \times K$ matrix inversion $Q_{t+1|t}^{-1}$ at each step, which collapses to $k \times k$ under a Kronecker covariance structure (see Time-Varying Parameter VAR). See State-Space Representation for the underlying Kalman filter.

Hit-and-Run Algorithm for Non-Conjugate Covariance Priors (Sun-Ni 2005)

When the prior on $\Sigma$ is the Yang-Berger (1994) reference prior $\pi_R(\Sigma) \propto |\Sigma|^{-1}\prod_{i<j}(d_i - d_j)^{-1}$ rather than the standard Jeffreys prior $|\Sigma|^{-(p+1)/2}$ , the conditional posterior $\pi(\Sigma|\phi, Y)$ has no inverse-Wishart form. Sun and Ni (2005) sample $\Sigma$ within the Gibbs sampler using a Metropolis hit-and-run step in the log-eigenvalue space:

Log-transform the current draw: decompose $\Sigma_{k-1} = OD O'$ (eigen) and form $\Sigma^* = O D^* O'$ where $d^*_i = \log(d_i)$ .
Select a random symmetric $p \times p$ direction matrix $V$ by drawing elements $z_{ij} \sim N(0,1)$ and normalizing: $v_{ij} = z_{ij}/(\sum_{l \le m} z_{lm}^2)^{1/2}$ .
Propose a new log-covariance $\Psi = \Sigma^*_{k-1} + \lambda V$ , $\lambda \sim N(0,1)$ ; decompose $\Psi = QC^*Q'$ .
Compute the log Metropolis ratio: $\alpha_k = \frac{T}{2}\sum_i(d^*_i - c^*_i) + \frac{1}{2}\text{tr}\!\left[(\exp\Sigma^*_{k-1})^{-1} - (\exp\Psi)^{-1}\right]S_k + \sum_{i<j}\!\left[\log(d^*_i - d^*_j) - \log(c^*_i - c^*_j)\right]$
Accept $\Sigma^* \leftarrow \Psi$ with probability $\min(1, e^{\alpha_k})$ ; recover $\Sigma_k = Q\,\text{diag}(e^{c^*_1},\ldots,e^{c^*_p})Q'$ .

The log-transform ensures $\Sigma$ remains positive definite throughout. Acceptance rate is ~32% in practice. The algorithm is nested as one block within a larger Gibbs loop over $(\phi, \Sigma)$ . See Noninformative Prior for VAR.

Griddy-Gibbs Sampler

When a full conditional $f(x_i \mid x_{-i})$ has no standard analytical form, the Griddy-Gibbs algorithm (Gelfand and Smith 1990) proceeds as follows:

Evaluate $f(x_i \mid x_{-i}^{(j)})$ over a grid of $G$ points $\{x_i^{(1)}, \ldots, x_i^{(G)}\}$ spanning the support.
Numerically integrate (e.g., trapezoidal rule) to compute the cumulative distribution function (CDF) $F(x_i \mid x_{-i}^{(j)})$ over the grid.
Draw $u \sim \mathrm{Uniform}(0,1)$ and invert the CDF by linear interpolation: $x_i^{(j)} = F^{-1}(u \mid x_{-i}^{(j)})$ .

The grid must cover the region of non-negligible probability; $G \approx 40$ is typically sufficient for unimodal conditionals. The key limitation is slow convergence if the conditional is multimodal — runs with different starting values help diagnose this.

Adaptive Rejection Sampling for Conjugate Hierarchical Hyperparameters (George-Makov-Smith 1993)

Griddy-Gibbs evaluates the log-conditional density on a numerical grid. When the conditional is known to be log-concave, Gilks-Wild (1992) adaptive rejection sampling (ARS) offers an exact alternative: it constructs a piecewise-exponential envelope using function evaluations and their derivatives, samples from the envelope by inversion, and accepts/rejects — adapting the envelope at each rejection, so the acceptance probability rises to 1 without tuning.

The hyperparameter bottleneck in conjugate hierarchical models. Consider an exponential-family hierarchy:

$\text{Data layer:} \quad x_i|\theta_i \sim f(x|\theta_i) \quad \text{(exponential family)}$ $\text{Prior layer:} \quad \theta_i|\alpha,\beta \sim \pi(\theta|\alpha,\beta) \quad \text{(natural conjugate, parameterized by }(x_0,n_0)\text{)}$ $\text{Hyperprior:} \quad \pi(\alpha,\beta)$

The $\theta_i$ blocks draw conjugately. The hyperparameter block has full conditional:

$\pi(\alpha,\beta|\theta_1,\ldots,\theta_p) \propto L(\alpha,\beta|\theta_1,\ldots,\theta_p)\cdot\pi(\alpha,\beta)$

where $L$ is the conjugate likelihood distribution — the likelihood of $(\alpha,\beta)$ treating $\theta_1,\ldots,\theta_p$ as i.i.d. draws from the conjugate prior:

$L(x_0,n_0|\theta_1,\ldots,\theta_p) \propto \exp\!\left[x_0\sum_i\theta_i - n_0\sum_i\psi(\theta_i) - p\varphi(x_0,n_0)\right] \tag{GMS eq. 21}$

Theorem 3.3 (George-Makov-Smith 1993): $L$ is log-concave in $(x_0,n_0)$ ; proper marginally for all $p\geq 1$ ; proper jointly iff $p\geq 2$ .

Gamma-Poisson case (Corollary 3.5): For $\theta_i|\alpha,\beta\sim\text{Gamma}(\alpha,\beta)$ :

$L^G(\alpha,\beta|\theta_1,\ldots,\theta_p) \propto \left(\frac{\beta^\alpha}{\Gamma(\alpha)}\right)^p\!\left(\prod_i\theta_i\right)^{\alpha-1}\exp\!\left(-\beta\sum_i\theta_i\right)$

$\alpha|\beta,\theta_{1:p}$ : log-concave → Gilks-Wild ARS (univariate, one evaluation per candidate)
$\beta|\alpha,\theta_{1:p} \sim \text{Gamma}(p\alpha+a,\, \sum_i\theta_i+b)$ → standard conjugate draw

Complete exact Gibbs cycle (pump failure, $p=10$ ):

$\theta_i|\alpha,\beta,x_i \sim \text{Gamma}(\alpha+x_i,\, \beta+t_i)$ — closed-form conjugate
$\alpha|\beta,\theta_{1:10}$ — ARS from $L^G(\alpha,\cdot)$
$\beta|\alpha,\theta_{1:10} \sim \text{Gamma}$ — closed-form conjugate

Every step is exact. No grid approximation or Metropolis acceptance/rejection step required. This is more efficient than Griddy-Gibbs when log-concavity holds, since ARS needs only the log-density and its gradient rather than a CDF computed on a fixed grid.

Contrast with Griddy-Gibbs: Griddy-Gibbs applies universally (any proper conditional) but approximates via CDF inversion on a finite grid. ARS applies only to log-concave conditionals but is exact and self-adapting. For the conjugate hierarchical hyperparameter problem, ARS is the preferred method when the conditional is one of the standard families (Gamma, Beta).

Rao-Blackwell Variance Reduction

Replacing the direct estimator of $E(X_i)$ : $\hat\mu_i = \frac{1}{k}\sum_{j=1}^k x_i^{(j)}$

with the Rao-Blackwellized estimator: $\hat\mu_i^c = \frac{1}{k}\sum_{j=1}^k E(X_i \mid x_{-i}^{(j)}) \tag{RB}$

always reduces variance: $\mathrm{Var}(\hat\mu_i^c) \leq \mathrm{Var}(\hat\mu_i)$ . With the Griddy-Gibbs sampler, $E(X_i \mid x_{-i}^{(j)})$ is the conditional mean over the grid — computed at each step anyway, so the RB estimator comes nearly for free.

The same principle extends to second moments: the within-group conditional variance $\mathrm{Var}(x_i \mid x_{-i})$ (known analytically from the conditional pdf) is averaged alongside the conditional mean to obtain a lower-variance estimate of $\mathrm{Var}(X_i)$ .

The principle extends equally to density estimation: instead of histogramming the sampled values, estimate $f(x)$ by averaging the conditional densities at the simulated companion values: $\hat{f}(x) = \frac{1}{m}\sum_{j=1}^m f(x \mid y_j) \tag{RB-density}$ Since $E[f(x|Y)] = \int f(x|y)f(y)\,dy = f(x)$ , equation (RB-density) is unbiased whenever $y_1,\ldots,y_m \sim f(y)$ . Casella and George (1992, eq. 2.9) demonstrate it on both continuous and discrete examples; the analytic justification is provided by Liu, Wong, and Kong (1994).

Extension to Accept-Reject and Metropolis (Casella-Robert 1996). The same ancillary-variable argument extends beyond Gibbs. In the Accept-Reject algorithm, the uniform variables $U_i \sim \text{Unif}[0,1]$ that decide whether to accept candidate $Y_i$ are ancillary. Rather than discarding rejected draws, one can compute Rao-Blackwellised weights: $p_i = P(Y_i \text{ accepted} \mid N, Y_1,\ldots,Y_N)$ via an $O(n^2)$ recurrence on $w_i = f(Y_i)/Mg(Y_i)$ , yielding estimator $\hat{\mu}^{RB} = \sum_i p_i h(Y_i)$ . mean squared error (MSE) reduction: 15–20% at acceptance rate 0.9, up to 60% at rate 0.3 — permanently, regardless of sample size. For the Metropolis algorithm, the analogous weights $\vartheta_i$ = expected number of times $Y_i$ appears in the chain (computed via a similar recursion) yield ~40–50% MSE reduction in the independent case and 7–25% in the random-walk case. A key empirical finding: a Rao-Blackwellised importance sampling estimator built on the Metropolis chain achieves 70–95% MSE reduction — substantially better than RB Metropolis itself — suggesting that the Metropolis mean is a poor estimation device (though the algorithm remains valuable as a sampler).

Blocked Gibbs Sampler

When two or more parameters $(\theta_i, \theta_j)$ have high posterior correlation, their marginal chains exhibit high autocorrelation and the effective sample size per draw drops sharply. The solution is to draw the correlated block jointly from the multivariate conditional $f(\theta_i, \theta_j \mid \theta_{-ij}, Y)$ :

Marginalize: compute $f(\theta_i \mid \theta_{-ij}, Y) = \int f(\theta_i, \theta_j \mid \theta_{-ij}, Y)\, d\theta_j$ numerically.
Draw $\theta_i$ from this marginal by Griddy-Gibbs.
Draw $\theta_j$ from the (now one-dimensional) full conditional $f(\theta_j \mid \theta_i, \theta_{-ij}, Y)$ .

Bauwens-Giot (1997) report 5–7× reduction in estimator variance for the price elasticity parameters ( $\delta_f$ , $\delta_d$ , posterior correlation $\approx -0.8$ ), with comparable CPU cost per unit of accuracy. The general principle is formalized in Liu, Wong, and Kong (1994).

Chib-Carlin (1999) — Blocked Gibbs for hierarchical longitudinal models. The Gaussian linear mixed model $y_i = X_i\beta + W_i b_i + \varepsilon_i$ has the standard 4-block Gibbs $(\beta, \{b_i\}, D^{-1}, \sigma^2)$ with severe mixing when variance components $D$ and $\sigma^2$ are both vague. Three progressively blocked algorithms address this: (1) Algorithm 2 — draw $\beta$ marginalized over $\{b_i\}$ using $y_i|\beta,\sigma^2,D\sim\mathcal{N}(X_i\beta, V_i)$ , $V_i=\sigma^2I+W_iDW_i'$ ; produces essentially iid $\beta$ draws (autocorrelation time $\kappa$ drops from 20 to 1 in the HIV CD4 example). (2) Algorithm 3 — evaluate the integrated likelihood $f(y|\sigma^2,D)$ via the Chib (1995) identity as a ratio of three Normal densities evaluated at the posterior mean $\beta^*$ : $f(y|\sigma^2,D)=\phi(\beta^*|\beta_0,B_0)\prod_i\phi(y_i|X_i\beta^*,V_i)/\phi(\beta^*|\hat\beta,B)$ ; update $(\sigma^2,D)$ jointly by M-H with multivariate- $t$ proposal; reduces $\kappa(D_{21})$ from 60 to 11. For binary longitudinal probit, marginalizing $Z_i$ over $b_i$ connects to the Chib-Greenberg (1998) multivariate probit; Algorithm 7 fully marginalizes over $b_i$ and $Z$ via the GHK likelihood to achieve near-iid $\beta$ draws at lag 10 (0.009 vs 0.934). Universal rule: always marginalize fixed effects over random effects before sampling. See Chib-Carlin (1999) and Bayesian Hierarchical Model.

Marginal Gibbs Sampler for Cointegrated VARs (Warne 2006)

When the VECM has a Kronecker-structured prior — $\Gamma_i|\Omega \sim \mathcal{N}(0,\Sigma_{\Gamma_i}\otimes\Omega)$ and $\mathrm{vec}(\alpha)|\beta,\Omega \sim \mathcal{N}(0,(\beta'\beta)^{-1}\otimes\lambda_\alpha^2\Omega)$ — the deterministic blocks $(\Phi,\Gamma)$ can be integrated out of the joint posterior analytically. This yields a marginal Gibbs sampler (Warne 2006, Proposition 5) that cycles over only three blocks:

$\Omega | \alpha, \psi, Y \sim \mathrm{IW}$ — the integration reduces the effective dimension of the Wishart posterior by $d + n(k-1)$ degrees of freedom (where $d$ = columns of $D_t$ , $n$ = VAR dimension, $k$ = lag order), producing a tighter posterior.
$\alpha | \psi, \Omega, Y \sim \mathcal{N}$ — closed-form conditional regression of loading matrix on cointegrating space.
$\psi | \alpha, \Omega, Y \sim$ matrix- $t$ — the Grassmann manifold parametrization $\Psi$ of the cointegrating space retains a matrix- $t$ full conditional (see Cointegration).

The cost is a $d + n(k-1)$ reduction in the Wishart degrees of freedom — an increase in posterior precision — because the integration over $(\Phi,\Gamma)$ concentrates mass. The full six-block sampler (Proposition 1) and the marginal three-block sampler yield the same marginal posterior for $(\alpha, \psi, \Omega)$ but the marginal sampler has lower autocorrelation per draw.

Rao-Blackwellization for the marginal likelihood. The marginal Gibbs output $\{\Psi^{(i)}\}$ is used to estimate $p(\tilde\alpha|Y,r)$ by averaging analytic Gaussian conditionals: $\hat{p}(\tilde\alpha|Y,r) = G^{-1}\sum_{i=1}^G p(\tilde\alpha|\Psi^{(i)},Y,r)$ . This is a Rao-Blackwell estimator for the posterior density at a point — the same principle as in the Bauwens-Giot and Savage-Dickey cases but applied to the Chib (1995) MLI denominator (see Marginal Data Density).

Extracting Samples from the Gibbs Sequence (Casella-George 1992)

Three main strategies for exploiting the Gibbs sequence, each with different trade-offs:

Final- $k$ method (Gelfand and Smith 1990): Run $m$ independent chains of length $k$ ; use the terminal value $X_k$ from each. Yields an approximately iid sample from $f(x)$ when $k$ is large enough. Diagnostic: compare density estimates from the $m$ chains and stop at the first $k$ at which they appear identical ("felt-tip pen test"). Simple to parallelize.
Thinning by $r$ (Geyer): Run one long chain; extract every $r$ th observation to reduce serial dependence. Advantage: less sensitivity to initial values. Risk: the chain may remain in a small subset of the sample space for an extended period (Gelman and Rubin 1991), which thinning does not cure.
All post-burn-in (George and McCulloch 1993): Use every $X_j$ for $j \geq k$ . Least wasteful for computationally expensive models; the empirical distribution still converges to $f(x)$ by ergodicity. The mixing rate of the chain determines the effective sample size per draw. In the SSVS variable-selection setting, all post-burn-in $\gamma$ draws are tabulated to estimate marginal inclusion probabilities $P(\gamma_i=1|Y)$ .

The Rao-Blackwell density estimator (RB-density) is applicable under any of these strategies and consistently outperforms direct histogram methods at no additional computational cost.

Subsampling is Never Beneficial (MacEachern-Berliner 1994)

Thinning — keeping only every $k$ th draw from a length- $nk$ run — never reduces the variance of posterior mean estimates and generally increases it. MacEachern-Berliner (1994) prove this for any stationary Markov chain, requiring no reversibility assumption and thus covering the standard fixed-scan Gibbs sampler.

The inefficiency proof. Partition the run into $k$ systematic subsample means: $\bar{x}_i = \frac{1}{n}\sum_{j=0}^{n-1} x_{i+jk}, \quad i = 1, \ldots, k$ By stationarity, all $k$ subsample means have the same variance: $\text{var}(\bar{x}_i) = \text{var}(\bar{x}_1)$ . The full-chain mean $\bar{x} = k^{-1}\sum_i \bar{x}_i$ satisfies: $\text{var}(\bar{x}) = \frac{\text{var}(\bar{x}_1)}{k} + \frac{1}{k^2}\sum_{i \neq j}\text{cov}(\bar{x}_i, \bar{x}_j) \leq \text{var}(\bar{x}_1)$ because $|\text{cov}(\bar{x}_i, \bar{x}_j)| \leq \text{var}(\bar{x}_1)$ . Equality holds only when the distribution is degenerate (constant). Hence the full chain weakly dominates any systematic subsample.

The stationary time series trap. The spectral variance formula $\text{var}(\bar{x}) = (\sigma^2/nk)\{1+2\sum(1-i/nk)\rho_i\}$ can suggest that thinning wins when $\rho_1 > 0.5$ : for an AR(1)-type autocorrelation sequence with $\rho_1 = \rho$ , $\rho_i = 0$ for $i > 1$ , and $k = 2$ , the thinned variance equals $\sigma^2/n$ while the full-chain variance equals $(\sigma^2/2n)(1 + 2\rho)$ , which exceeds $\sigma^2/n$ when $\rho > 0.5$ . But stationarity with this autocorrelation structure requires $|\rho| \leq 0.5$ — the scenario never arises. The constraints on the autocorrelation function of any stationary process close off every parameter combination that would make thinning beneficial.

Extension beyond Geyer (1992). Geyer's (1992) analogous result required chain reversibility (detailed balance). MacEachern-Berliner require only stationarity, since Tierney (1992) noted that the standard fixed-scan Gibbs sampler is typically non-reversible.

Preferred alternatives. Batching (block-mean variance estimator) and spectral variance estimation (see Spectral Variance Estimator above) use the full chain and are uniformly preferred for variance estimation. The only legitimate reason to thin is computational: limited storage capacity before processing, or when processing cost per draw exceeds generation cost — a scenario more relevant in 1994 than today.

Predictive Distributions via MCMC (Scollnik ~1996)

A primary actuarial output is the posterior predictive distribution for future data $Y_2$ given observed data $Y_1$ :

$f(Y_2 \mid Y_1) = \int f(Y_2 \mid Y_1, \theta)\, f(\theta \mid Y_1)\, d\theta$

The MCMC estimate follows directly from the ergodic theorem (eq. 2.2): average $f(Y_2 \mid Y_1, \theta^{(i)})$ over posterior draws $\theta^{(i)}$ . This is the same principle as the Rao-Blackwell density estimator (RB-density) — conditional densities averaged over sampled parameters — and requires no extra computation beyond the Gibbs run itself.

Posterior inference on any transformation $w = w(\theta)$ is equally free: compute $w^{(i)} = w(\theta^{(i)})$ and form the empirical distribution.

Data Augmentation for Grouped Data (Scollnik ~1996, Example 3)

When observations are grouped (only class frequencies observed, not exact values), a data augmentation approach adds the unobserved exact values as latent auxiliary variables, making the full conditionals tractable. For $n$ losses from Pareto( $\theta, \lambda$ ) observed in $C$ class intervals with counts $f_c$ :

Augment: introduce latent exact losses $x_i$ for each observation; their conditional distribution given $\theta$ and the class interval is a truncated Pareto on that interval.
Update $\theta$ : with exact values in hand, the conjugate gamma posterior is immediate — $\theta \mid x, \lambda \sim \text{Gamma}(n + \alpha,\; \beta - n\ln\lambda + \sum_i \ln(\lambda + x_i))$ .
Iterate: the Gibbs chain over $(\theta, x_1, \ldots, x_n)$ converges quickly (empirically, ~10 iterations suffice for 25 observations in 4 classes).

Side inferences on tail events — e.g., P(at least one open-interval loss $> L$ ) — are read directly from the simulated augmented data, requiring no additional computation.

The augmentation idea generalizes: any latent-variable model (censored data, mixture models, state-space models) can be handled by adding the missing data as a Gibbs block.

Data Augmentation for Variable Selection (George-McCulloch 1993)

Stochastic Search Variable Selection (SSVS) casts subset selection as a Gibbs problem over $(\beta, \sigma^2, \gamma)$ , where $\gamma \in \{0,1\}^p$ is a latent binary inclusion vector added as an augmentation block. See Variable Selection for the full model; the Gibbs structure is:

Block 1 — $\beta^j$ : Conjugate Normal draw given current $(\sigma^{j-1}, \gamma^{j-1})$ — the mixture prior $\beta|\gamma \sim \mathcal{N}_p(0, D_\gamma R D_\gamma)$ yields a closed-form posterior: $\beta^j \sim \mathcal{N}_p\!\left(A_\gamma(\sigma^{j-1})^{-2}X'X\hat\beta_{LS},\; A_\gamma\right), \quad A_\gamma = \left((\sigma^{j-1})^{-2}X'X + D_\gamma^{-1}R^{-1}D_\gamma^{-1}\right)^{-1}$

Block 2 — $\sigma^j$ : Conjugate Inverse-Gamma draw given current $(\beta^j, \gamma^{j-1})$ .

Block 3 — $\gamma_i^j$ : Independent Bernoulli draws. For each $i$ , the inclusion probability is $a/(a+b)$ where $a = f(\beta_i^j|\gamma_i=1)\cdot p_i$ and $b = f(\beta_i^j|\gamma_i=0)\cdot(1-p_i)$ — evaluating the spike and slab Normal density at the already drawn $\beta_i^j$ . Crucially, $Y$ enters only through $\beta^j$ , not directly.

The all-post-burn-in strategy reads off the empirical frequency of $\gamma_i=1$ across draws as an estimate of $P(\gamma_i=1|Y)$ , navigating $2^{60}$ model space without enumeration.

MN-GARCH Gibbs Sampler (Bauwens-Rombouts 2005)

The Bauwens-Rombouts (2005) sampler for the Mixture-of-Normals GARCH model illustrates two non-standard Gibbs patterns within a GARCH context:

Data augmentation for the mixture: Latent state variables $S_t \in \{1,\ldots,K\}$ are added as a Gibbs block. Conditional on states, each GARCH component's likelihood factors, converting the mixture problem into $K$ separate conditional problems. The state draw is multinomial: $P(S_t = k | \pi, \mu, \theta, y_t) \propto \pi_k \phi(y_t; \mu_k, h_{k,t})$ .

Griddy-Gibbs for GARCH parameters: The GARCH parameters $(\omega_k, \alpha_k, \beta_k)$ per component have no conjugate conditional because GARCH variances $h_{k,t}$ depend on past observations recursively. The griddy-Gibbs algorithm evaluates the log-conditional density on a 3-D grid, numerically integrates to form the CDF, and inverts via linear interpolation. Positivity constraints $\omega_k > 0$ , $\alpha_k \geq 0$ , $\beta_k \geq 0$ , and the stationarity constraint $\alpha_k + \beta_k < 1$ (for stationary components) bound the grid.

Dirichlet conjugate for mixing weights: $\pi | S_T \sim \text{Dirichlet}(a_{k0} + n_k)$ where $n_k = \#\{t: S_t = k\}$ . This is the same conjugate update as in the i.i.d. mixture case; conditioning on $S_T$ makes it valid even in the GARCH setting.

Constrained mean block: The $K-1$ free means $(\mu_1, \ldots, \mu_{K-1})$ are drawn jointly from a multivariate Normal that enforces the zero-mean constraint $\sum_k \pi_k \mu_k = 0$ via a $(K-1)\times(K-1)$ linear substitution.

See GARCH and BEKK-GARCH and Mixture of Normals for the model context.

Three-Block SUR Gibbs with Cross-Set Coefficient Pooling (Griffiths-Valenzuela 2002)

When $H$ groups of SUR equations share a common coefficient vector $\eta$ but have group-specific parameters $\Theta_h$ and covariance matrices $\Omega_h$ , the fully conjugate Gibbs sampler decomposes into three natural blocks:

Model: $Y_h = Z_h\Theta_h + X_h\eta + e_h$ , $h = 1,\ldots,H$ ; $e_h \sim \mathcal{N}(0, \Omega_h \otimes I_{M_h})$ ; each set $h$ has $M_h$ observations and $N_h$ equations.

Block 1 — $\Omega_h \mid \Theta_h, \eta, Y \sim \mathrm{IW}(A_h, M_h)$ : The residual within set $h$ after removing both set-specific and common effects is $\tilde{e}_h = \tilde{Y}_h - \tilde{Z}_h\Theta_h - \tilde{X}_h\eta$ , and $A_h = \tilde{e}_h'\tilde{e}_h$ (reshaped to $N_h \times N_h$ ). One Inverted-Wishart draw per set, all $H$ draws are independent given the other blocks.

Block 2 — $\Theta_h \mid \Omega_h, \eta, Y \sim \mathcal{N}(\hat\Theta_h, [Z_h'(\Omega_h^{-1} \otimes I_{M_h})Z_h]^{-1})$ : Within set $h$ , partial out the common $\eta$ contribution from the dependent variable, then run a standard GLS draw. The precision matrix $Z_h'(\Omega_h^{-1} \otimes I_{M_h})Z_h$ is block-diagonal within the set.

Block 3 — $\eta \mid \{\Omega_h\}, \{\Theta_h\}, Y \sim \mathcal{N}(\hat\eta, W^{-1})$ : The precision matrix pools information across all $H$ sets: $W = \sum_{h=1}^H X_h'(\Omega_h^{-1} \otimes I_{M_h})X_h, \qquad Q = \sum_{h=1}^H X_h'(\Omega_h^{-1} \otimes I_{M_h})(\tilde{y}_h - \tilde{Z}_h\Theta_h)$ $\hat\eta = W^{-1}Q$

This is the novel block: it is the GLS estimator for $\eta$ that treats each set as contributing observations to a pooled regression, weighted by its own $\Omega_h^{-1}$ . Compared to Percy (1992) and Chib-Greenberg (1995b), which handle single-set SUR, this block aggregates cross-set information in one conjugate step.

Key distinction from Baltagi (1995): Baltagi assumes the same $\Omega_h = \Omega$ for all sets; Griffiths-Valenzuela allow distinct unrestricted $\Omega_h$ per set, which is more flexible but requires separate Block 1 draws for each $h$ .

Applications (Griffiths-Valenzuela 2002):

ELES household equivalence scales: $H=8$ demographic types, $N_h=11$ commodities, 5,532 Australian households; 18,000 draws (3,000 burn-in)
Merino wool translog cost: $H=23$ years (sets), time-varying intercepts $\Theta_h$ , common $\alpha_i$ , $\alpha_{ij}$ technology parameters; 23,000 draws (3,000 burn-in)

See Equivalence Scale and Seemingly Unrelated Regression.

Sparse Factor Analytic Probit Gibbs Sampler (Hahn-Carvalho-Scott 2012)

Hahn, Carvalho, and Scott (2012) construct a 12-step conjugate Gibbs sampler for the sparse factor analytic probit model applied to high-dimensional binary outcomes. The sampler extends Chib-Greenberg (1998) by replacing the unstructured correlation matrix with a sparse factor decomposition.

Model: $z_i \sim \mathcal{N}(\alpha + \lambda\gamma_i + Bf_i, I_p)$ ; $y_{ij} = \mathbf{1}[z_{ij} > 0]$ ; $\Sigma = BVB^\top + I$ ; $B$ lower triangular ( $b_{jj} > 0$ , $b_{js} = 0$ for $s > j$ ); $f_i \sim \mathcal{N}(0, V)$ ; $V$ free $k \times k$ positive-definite.

Spike-and-slab prior on loadings: $(b_{js} \mid \nu_s, q_s) \sim q_s \cdot \mathcal{N}(0, \nu_s) + (1-q_s) \cdot \delta_0(b_{js})$ ; $\nu_s \sim \mathrm{IG}$ ; $q_s \sim \mathrm{Beta}(1,1)$ .

Selected Gibbs blocks:

Step	Variable	Conditional	Form
3	$\gamma_i$	$z_i, \alpha, \lambda, B, f_i, V$	Truncated Normal: $\mathbb{R}^+$ for Republicans, $\mathbb{R}^-$ for Democrats
5	$f_i$	$z_i, \alpha, \lambda, \gamma_i, B, V$	$\mathcal{N}(m, M)$ ; $M = (V^{-1} + B^\top B)^{-1}$ ; $m = M \cdot B^\top(z_i - \alpha - \lambda\gamma_i)$
6	$V$	$\{f_i\}, V_0$	$\mathrm{IW}(V_0 + FF^\top, n + h)$
7	$b_{js}$	$z, f, \alpha, \lambda, \gamma, \nu_s, q_s$	Spike-and-slab: compute Bayes factor $\hat{r}_{js}$ ; $\hat{q}_{js} = \hat{r}_{js}/[(1-q_s)/q_s + \hat{r}_{js}]$
10	$z_{ij}$	$y_{ij}, \alpha, \lambda, \gamma_i, B, f_i$	Truncated $\mathcal{N}$ : $(0,\infty)$ if $y_{ij}=1$ , $(-\infty,0)$ if $y_{ij}=0$

Key design choices:

Free $V$ vs. orthogonal factors. Setting $V = I_k$ forces all inter-factor correlation into $B$ , producing a denser loading matrix. A free $V$ absorbs inter-factor correlation, so $B$ requires fewer non-zero entries. The lower-triangular constraint (Geweke-Zhou 1996) prevents the rotational non-identifiability that would otherwise make $V$ and $B$ separately unidentified.
Multiplicity correction. The $\mathrm{Beta}(1,1)$ prior on $q_s$ (inclusion probability for factor column $s$ ) implements the Scott-Berger (2006, 2010) automatic multiplicity correction: when the number of candidate loadings $p$ is large, the marginal prior $p(\text{\# included}) \propto 1/(\text{number of models of that size})$ , which down-weights large models.

Simulation results (Table 1, n=50):

Model	p=20 Stein loss	p=100 Stein loss
Sparse-6-factor (k=3 true)	3.0	43.4
Factor (no sparsity)	4.3	136.7
Wishart	11.9	503.1

Senate application (1949–2009): Three-factor model on 30 closest roll-call votes per term. Partisanship factor $\gamma_i$ extracted as party-aligned dimension; dimensionality of bills stable across decades; partisanship low in late 1960s–early 1970s, increased sharply from 1979.

See Sparse Factor Probit, Multivariate Probit, and Variable Selection.

Data Augmentation for Scale Mixtures: Student-t Disturbances (Geweke 1993)

Geweke (1993) shows that a linear regression with i.i.d. Student-t( $\nu$ ) errors can be sampled by a four-block conjugate Gibbs sampler once the latent scale weights $\omega_i$ are included. The equivalence $y_i \mid x_i, \omega_i \sim \mathcal{N}(x_i'\beta, \sigma^2\omega_i)$ with $\nu/\omega_i \sim \chi^2(\nu)$ renders the following full conditionals all conjugate:

Block	Variable	Conditional	Form
1	$\omega_i$	$(\beta,\sigma,\nu)$	IG: $[(\sigma^{-2}u_i^2+\nu)/\omega_i]\sim\chi^2(\nu+1)$
2	$\sigma^2$	$(\beta,\omega)$	IG: $[\sum_i u_i^2/\omega_i]/\sigma^2\sim\chi^2(n)$
3	$\beta$	$(\sigma,\omega)$	Normal: GLS $\hat\beta(\omega)=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y$
4	$\nu$ (optional)	$(\beta,\sigma,\omega)$	Non-standard kernel; exponential prior recommended

Posterior odds ratios between competing degrees-of-freedom values are available as a by-product of the Gibbs output, as an expectation under the $\nu^{(2)}$ -posterior (eq. 16). Convergence proofs (Theorems 4–5) are provided for fixed $\nu$ and for the exponential-prior case.

Connection to SV models. The cross-sectional $\omega_i$ weights here are the static predecessor of the time-varying $h_t$ in the Jacquier-Polson-Rossi (1994) stochastic volatility model: JPR's data-augmentation with a latent variance path is the dynamic extension of Geweke's latent scale weights. See Geweke (1993), Mixture of Normals, and Stochastic Volatility.

Truncated Normal Gibbs for Inequality-Constrained Regression (Geweke 1996)

Geweke (1996) shows how to impose linear inequality constraints $a \leq D\beta \leq w$ on the coefficients of a normal linear regression entirely within a Gibbs sampler, without the accept-reject inefficiency of the crude frequency simulator.

Key reparameterization. Let $z = D(\beta - b)$ , $a^* = D(a-b)$ , $w^* = D(w-b)$ , $T = \sigma^2 D(X'X)^{-1}D'$ . The constraint region $a \leq D\beta \leq w$ becomes $a^* \leq z \leq w^*$ coordinate by coordinate. The two full conditionals are:

$\sigma \mid z$ : $(y - X\beta)'(y - X\beta)/\sigma^2 \sim \chi^2(T)$ — same chi-squared as the unconstrained model (eq. 3.2).
$z_j \mid \sigma, z_{i\neq j}$ : univariate Normal truncated to $[a_j^*, w_j^*]$ , with mean $b_j + \sum_{i\neq j} c_{ij}z_i$ and variance $h_j^2 = (T_{jj})^{-1}$ (eq. 3.3).

The critical insight is that rotating the problem by $D$ converts joint constraints into independent per-coordinate bounds, so each $z_j$ 's full conditional is a univariate (not multivariate) truncated normal — directly samplable via inverse-CDF. The chain converges in distribution to the constrained posterior because there is positive probability of moving to any positive-posterior region in one step.

Computational advantage. When the constraint probability $p_{2|1}$ is small (e.g., $9\times10^{-5}$ for the automobile sales example with k=11 sign constraints), the crude accept-reject simulator requires hours for 1,000 valid draws; the Gibbs sampler achieves equivalent accuracy in 49 minutes using every-100th draw.

Extensions. The same template applies to any model where a parameter block's conditional distribution is truncated multivariate normal: SUR models, censored regression, multinomial probit, and scale-mixture regression (Student-t with inequality constraints — the $z_j$ draws inherit the scale weights from Geweke 1993). See Geweke (1996) and Multinomial Probit.

Mixture of Normals Probit Gibbs Sampler (Geweke-Keane 1997)

Geweke and Keane (1999) construct a fully conjugate Gibbs sampler for binary probit with a $m$ -component Gaussian mixture disturbance. The augmented state at each iteration consists of: the Albert-Chib (1993b) latent utilities $\{\tilde{y}_t\}$ and a new set of component-assignment indicators $\{L_t\} \in \{1,\ldots,m\}^T$ . The six blocks, all with standard-form conditionals, are:

Block	Variable	Conditional	Form
1	$\tilde{y}_t$	$d_t, L_t, \beta, \alpha, h$	Truncated Normal on $(0,\infty)$ or $(-\infty,0)$ ; mean $\alpha_{L_t}+\beta'x_t$ , precision $h_{L_t}$
2	$L_t$	$\tilde{y}_t, \beta, \alpha, h, p$	Multinomial: $P(L_t=j)\propto p_j\exp[-0.5h_j(\tilde{y}_t-\alpha_j-\beta'x_t)^2]$
3	$\beta$	$\tilde{y}, L, \alpha, h$	Normal: precision-weighted regression with $\sum_t h_{L_t} x_t x_t'$
4	$\alpha$	$\tilde{y}, L, \beta, h$	Normal with ordering constraint; draws via Geweke (1991)/(1995) truncated-Normal algorithm
5	$h_j$ ( $j\neq j^*$ )	$\tilde{y}, L, \beta, \alpha$	$s_j^2 h_j \sim \chi^2(\nu_j+T_j)$ ; $s_j^2$ updated by component SSR
6	$p$	$L$	Dirichlet with counts $r_j + T_j$

Dual augmentation principle. Block 1 is the Albert-Chib (1993b) latent-utility device; block 2 is the Chib (1996)-style component-assignment device. Together they handle both the binary truncation and the mixture non-Gaussianity within a single unified Gibbs cycle. The same pattern (latent utility + component label) generalizes to multinomial probit with mixture disturbances (Geweke-Keane-Runkle 1994, 1997).

Marginal likelihood. The Gelfand-Dey (1994) estimator, using a truncated multivariate Normal function $f(\cdot)$ centered at posterior mean of $(\beta,\alpha,h,p)$ , avoids the dimension explosion that would arise from including $(\tilde{y},L)$ in the parameter space; the analytically marginalized likelihood (eq. 2.15) is $p(d|\beta,\alpha,h,p,X) = \prod_t [d_t\sum_j p_j\Phi(h_j^{1/2}(\alpha_j+\beta'x_t)) + (1-d_t)(1-\sum_j p_j\Phi(\ldots))]$ .

See Geweke-Keane (1999), Binary Probit, and Mixture of Normals.

Data Augmentation for Panel Data and Mixture Shocks (Geweke-Keane 2000)

Geweke and Keane (2000) illustrate two non-standard Gibbs patterns in a dynamic panel earnings model fitted to the PSID:

Unbalanced panel augmentation. Individuals enter and exit the panel at different ages, breaking the balanced-panel structure required for tractable AR-based conditionals. Unobserved period earnings $y_{it}$ are added as latent variables; their full conditional is Gaussian from the model's AR structure, enabling a complete balanced panel at each Gibbs iteration.

Mixture-of-normals block. First-period shocks $e_{i1}$ and transitory increments $g_{it}$ each follow a 3-component mixture $\sum_{k=1}^3 p_k\, \mathcal{N}(a_k,\, h_k^{-1})$ . Latent component indicators $s_{it} \in \{1,2,3\}$ are added as a Gibbs block; conditional on $s_{it}$ , all parameters have conjugate Gaussian or Gamma conditionals. Proper Gamma priors on precisions $h_k$ are essential: the mixture likelihood is unbounded as $h_k \to \infty$ for any non-empty component unless regularized. The normalization $a_2 = 0$ , $a_1 \leq a_2 \leq a_3$ prevents label switching.

The 8-block earnings sampler (plus marital status and individual effects $q_i$ ) achieves RNE of 0.05–1.0; the full PSID sample requires ~332 sec per iteration and 2,500 sweeps. See Mixture of Normals for the general framework.

Auxiliary Variable Gibbs Sampler for Non-Conjugate Models (Damien-Wakefield-Walker 1999)

When a full conditional takes the non-standard form $f(x) \propto \pi(x)\prod_{i=1}^N l_i(x)$ — where $\pi$ is a proper density and each $l_i$ is non-negative and invertible — Damien, Wakefield, and Walker (1999) show how to construct a Gibbs sampler with only uniform and truncated-standard conditionals.

Theorem 1. Introduce latent variables $U_i \sim \mathrm{Uniform}(0, l_i(x))$ for $i=1,\ldots,N$ . The joint density of $(X, U_1,\ldots,U_N)$ is: $f(x, u_1,\ldots,u_N) \propto \pi(x)\prod_{i=1}^N \mathbf{1}[0 < u_i < l_i(x)]$

The full conditionals are:

$U_i \mid x \sim \mathrm{Uniform}(0,\, l_i(x))$ — trivially sampled
$X \mid u_1,\ldots,u_N \sim \pi(x)$ restricted to $A_u = \{x : l_i(x) > u_i,\; i=1,\ldots,N\}$ — a draw from $\pi$ truncated to an intersection of intervals

The key requirement is that the invertibility of each $l_i$ gives an analytical expression for $A_u$ , making the truncated draw tractable.

"Always moves" property. The auxiliary variable sampler is equivalent to a Metropolis-Hastings independence chain with proposal $\pi(\cdot)$ but one that always moves — the chain never rejects. This is because sampling $U$ before $X$ conditions the proposal on falling in the acceptance region of the implied MH step. This dominates the corresponding MH independence chain for any estimand (Peskun 1973 optimality argument).

Choosing latent variables. The decomposition $f = \pi \cdot \prod l_i$ is not unique; this flexibility is exploited to obtain computable truncation sets. For example:

$l(x) = \exp\{-\exp(x)\}$ : introduce $U$ such that $f(x,u) \propto \exp(-u)\mathbf{1}\{u > \exp(x)\}\pi(x)$ , giving $f(U|x) \sim \mathrm{Exp}(1)$ restricted to $(\exp(x), \infty)$ and $f(X|u) \sim \pi$ truncated to $\{x: x < \log(u)\}$ .
Logistic likelihood: introduce one $U_i$ and one $V_i$ per observation; truncation set becomes an interval in $x$ defined by min/max over $S = \{i: y_i=1\}$ and $R = \{i: y_i=0\}$ .

Applications in the paper:

Poisson-log-normal model: $f(x|r) \propto \exp\{rx - \exp(x) - 0.5x^2\}$ → latent exponential $U$ → truncated Normal for $X$ .
Bernoulli-logistic regression: two sets of latents $U_i, V_i$ per observation → truncated Normal for $\beta$ on a computable interval $A_w$ .
Weibull proportional hazards: exponential latents for the survival term → tractable interval bounds for both $\beta$ and the baseline hazard $a$ .
Generalized linear mixed models (GLMMs): two-latent-per-observation scheme for both the regression block $\beta$ and random effects $b_i$ ; all full conditionals become truncated Normal.
Non-linear mixed models (logistic growth): dimension-by-dimension sampling from $\pi(\theta_k|\theta_{-k})$ restricted to $A^k_u$ when the multivariate set $A_u$ is intractable.

Numerical comparison (Lindstrom-Bates 1990 orange-tree growth data): 7 trees, 5 observations each, 3-dimensional random effect $\theta_i$ . Auxiliary variable Gibbs vs. random-walk MH for 10,000 iterations:

Lag-1 autocorrelation: 0.75 (auxiliary) vs. 0.79 (MH)
Autocorrelation at lag 40: ≈0 (auxiliary) vs. 0.39 (MH)
Run times approximately equal; MH required tuning of the proposal scale.

Caveats. The decomposition $f = \pi \cdot \prod l_i$ must be discovered case-by-case; there is no automatic algorithm. When an efficient rejection sampler or MH proposal already exists, the auxiliary variable approach may not improve efficiency. The "always moves" property does not guarantee rapid mixing if the truncation set $A_u$ is small relative to the support of $\pi$ .

See Damien-Wakefield-Walker (1999), Paul Damien, Jon Wakefield, Stephen Walker, and Metropolis-Hastings Algorithm.

Convergence Theory (Tierney 1994)

Tierney (1994) establishes formal convergence properties using general state space Markov chain theory (Nummelin 1984). The key hierarchy:

Theorem 1 (Ergodic convergence). If the chain has invariant distribution $\pi$ , is $\pi$ -irreducible, and is aperiodic, then $\pi$ is the unique invariant distribution and for $\pi$ -almost all starting points $x$ : $\|P^n(x,\cdot) - \pi(\cdot)\| \to 0$

Geometric ergodicity via drift condition. Suppose there exist a function $V: E \to [1,\infty)$ , a "small" set $C$ (where $P(x,\cdot) \geq \epsilon\nu(\cdot)$ uniformly), and constants $\lambda < 1$ and $b < \infty$ such that: $PV(x) \leq \lambda V(x) + b\cdot\mathbf{1}_C(x) \quad \forall x$

Then the chain is geometrically ergodic: $\|P^n(x,\cdot) - \pi\| \leq M(x)\rho^n$ for some $\rho < 1$ (Proposition 1 / Tierney's formalization of Nummelin). This is the standard sufficient condition for central limit theorem (CLT) validity.

Theorem 4 (CLT for geometric ergodic chains). If the chain is geometrically ergodic and $g \in L^{2+\varepsilon}(\pi)$ , then: $\sqrt{n}\left(n^{-1}\sum_{j=1}^n g(X_j) - E_\pi[g]\right) \xrightarrow{d} \mathcal{N}(0,\sigma^2)$

where the asymptotic variance is the spectral density at zero: $\sigma^2 = \mathrm{Var}_\pi(g) + 2\sum_{j=1}^\infty \mathrm{Cov}_\pi(g(X_0),\, g(X_j))$

Geometric ergodicity guarantees this series converges absolutely. For uniform ergodicity (stronger: $\|P(x,\cdot)-\pi\|$ bounded uniformly), $L^2(\pi)$ functions suffice (Theorem 5).

Practical implication. For a Gibbs or MH sampler to yield valid confidence intervals (via batch means or spectral variance estimation), geometric ergodicity must hold. For posteriors with sub-exponential tails, the drift condition is typically satisfied. For posteriors with heavy tails or bounded supports, uniform ergodicity often holds directly.

Peskun (1973) optimality of MH acceptance ratio. Among all Metropolis-type acceptance functions $\alpha(x,y)$ that leave $\pi$ invariant, the Hastings form $\alpha(x,y) = \min\{\pi(y)q(y,x)/[\pi(x)q(x,y)],\, 1\}$ minimises the rejection probability for any pair $(x,y)$ — and therefore minimises the asymptotic variance $\sigma^2$ for any estimand. Barker's (1965) acceptance form $\pi(y)q(y,x)/[\pi(x)q(x,y) + \pi(y)q(y,x)]$ is strictly worse by this criterion.

Metropolis-Hastings Algorithm (Jackman 2000)

The Metropolis-Hastings (MH) algorithm generalizes Gibbs sampling to cases where a full conditional cannot be sampled directly. At each iteration, a candidate $\theta^*$ is drawn from a proposal density $q(\theta^*|\theta^{t-1})$ , then accepted with probability:

$r = \min\!\left(1,\; \frac{p(\theta^* \mid Y)\, q(\theta^{t-1} \mid \theta^*)}{p(\theta^{t-1} \mid Y)\, q(\theta^* \mid \theta^{t-1})}\right)$

If rejected, the chain stays at $\theta^{t-1}$ . The Gibbs sampler is a special case: when $q(\theta^*|\theta^{t-1}) = p(\theta^* \mid \theta^{t-1}_{-i}, Y)$ (a full conditional), the ratio $r$ equals 1 everywhere and every candidate is accepted.

Practical choices for q:

Random-walk proposal: $\theta^* = \theta^{t-1} + \varepsilon$ , $\varepsilon \sim \mathcal{N}(0, \tau^2)$ . Acceptance rate $\approx 0.44$ for scalar $\theta$ , $\approx 0.23$ for high-dimensional $\theta$ .
Independence sampler: $q(\theta^*|\theta^{t-1}) = q(\theta^*)$ — draws are independent of the current state. Efficient when $q$ approximates $p(\theta|Y)$ well; degenerates otherwise.

When MH steps are nested inside a larger Gibbs cycle, the resulting algorithm is called Metropolis-within-Gibbs (or hybrid MH-Gibbs). See Jackman (2000) for worked examples in probit and AR(1) models.

Convergence Diagnostics

Gelman-Rubin Potential Scale Reduction (Gelman and Rubin 1992)

Run $m \geq 2$ independent chains of length $2n$ from widely dispersed starting values; discard the first $n$ draws as burn-in. For a scalar parameter $\psi$ :

$B/n$ = between-chain variance of chain means
$W$ = average within-chain variance

The potential scale reduction factor $\hat{R} = \sqrt{\hat{V}/W}$ where $\hat{V} = (n-1)W/n + B/n$ . At convergence, $\hat{R} \to 1$ ; values substantially above 1 (say $> 1.1$ ) indicate insufficient mixing. The intuition: if all chains have converged to the same distribution, the between-chain variance should equal the within-chain variance.

Raftery-Lewis Chain Length Determination (Raftery and Lewis 1992)

Given a target quantile $q$ , tolerance $r$ , and confidence level $s$ , the method computes the required burn-in $M$ , thinning interval $k$ , and total chain length $N$ from a short pilot run (~2,000 draws).

Step 1 — Binarize. Set $Z_j = \mathbf{1}[\theta^{(j)} \leq \hat{t}]$ where $\hat{t}$ is the pilot-estimated $q$ -th quantile. Under stationarity, $\{Z_j\}$ is approximately a first-order two-state Markov chain.

Step 2 — Fit transitions. Estimate $\hat\alpha = n_{01}/(n_{00}+n_{01})$ and $\hat\beta = n_{10}/(n_{10}+n_{11})$ , where $n_{ij}$ counts $i \to j$ transitions. Check first-order adequacy via a $G^2$ deviance test vs. a second-order model; if rejected, thin by $k$ (use only every $k$ -th draw) until first-order holds.

Step 3 — Required burn-in. $M \approx \left\lceil \frac{\log\!\bigl(\varepsilon(\hat\alpha+\hat\beta)/\max(\hat\alpha,\hat\beta)\bigr)}{\log|1-\hat\alpha-\hat\beta|} \right\rceil$ where $\varepsilon$ is a small constant (e.g. 0.001) controlling closeness to stationarity.

Step 4 — Required chain length. $N_{\min} \approx k \cdot \frac{2\, z_{(1+s)/2}^2\, (\hat\alpha+\hat\beta-\hat\alpha^2-\hat\beta^2)}{r^2\, q(1-q)\, (\hat\alpha+\hat\beta)^2}$

For low-autocorrelation chains, $N \approx 3{,}000$ – $5{,}000$ is typical. Slowly mixing chains can require $N > 50{,}000$ .

More targeted than Gelman-Rubin: rather than assessing whether convergence has been reached, it tells you how many iterations you need for a given inferential precision. The method is quantile-specific — apply to the most demanding quantile when multiple are of interest.

Zellner-Min (1995) Convergence Criteria: GSC²

Gelman-Rubin and Raftery-Lewis diagnostics detect whether the chain has mixed — they do not detect convergence to an incorrect distribution. If the GS is trapped near one mode of a multimodal posterior, both diagnostics can indicate convergence while the estimated marginals are wrong.

Zellner and Min (1995) exploit the fact that when the joint posterior $p(\alpha,\beta \mid D) \propto \pi(\alpha,\beta \mid I_0)\, l(\alpha,\beta \mid y)$ is known analytically up to a constant, two factorizations must hold simultaneously:

$p(\alpha \mid D)\,p(\beta \mid \alpha, D) = p(\beta \mid D)\,p(\alpha \mid \beta, D) \tag{ZM-identity}$

This identity can be evaluated pointwise using GS-estimated marginals $\hat{p}_N(\alpha_i \mid D)$ , $\hat{p}_N(\beta_i \mid D)$ alongside the known conditional densities $p(\beta \mid \alpha, D)$ and $p(\alpha \mid \beta, D)$ .

Three GSC² criteria evaluated at parameter-space points $(\alpha_i, \beta_i)$ :

Criterion	Formula	Converged value
ARC² (Anchored Ratio)	$\hat\theta_{12}(N) = \hat{p}_N(\alpha_1\\|D)\,p(\beta_1\\|\alpha_1,D) / [\hat{p}_N(\alpha_2\\|D)\,p(\beta_2\\|\alpha_2,D)]$ vs. exact $\theta_{12}$	$\hat\theta_{12} \approx \theta_{12}$
DC² (Difference)	$\hat\eta_i(N) = \hat{p}_N(\alpha_i\\|D)\,p(\beta_i\\|\alpha_i,D) - \hat{p}_N(\beta_i\\|D)\,p(\alpha_i\\|\beta_i,D)$	$\hat\eta_i \approx 0$
RC² (Ratio)	$\hat\gamma_i(N) = \hat{p}_N(\alpha_i\\|D)\,p(\beta_i\\|\alpha_i,D) \;/\; \hat{p}_N(\beta_i\\|D)\,p(\alpha_i\\|\beta_i,D)$	$\hat\gamma_i \approx 1$

Inference on $\hat\eta_i$ . Under large-sample theory, $\hat\eta_i(N) \approx \mathcal{N}(\eta_i,\, \hat\phi_i^2)$ where the variance $\hat\phi_i^2$ (equation 12) accounts for GS autocorrelation via the sample autocovariances of the GS draws. Posterior odds for $H_0: \eta_i = 0$ vs. $H_A: \eta_i \neq 0$ : $K_{OA} = \sqrt{\pi N/2}\,(1 + z_i^2/N)\exp\{-z_i^2/2\}, \qquad z_i = \hat\eta_i / \hat\phi_i$ This uses a Cauchy prior on $\eta_i$ under $H_A$ (Jeffreys 1967). For the SUR example, $K_{OA}$ rises from 10.98 at $N=100$ to 118.68 at $N=9{,}000$ , strongly favoring $H_0$ and confirming correct convergence.

Bimodal failure detection. For the mixture $p = 0.1\, p_A + 0.9\, p_B$ (modes at $(0,5)$ and $(0,-5)$ ), the GS gets stuck near $(0,5)$ . Although the chain appears stationary, the DC² criterion gives $\hat\eta = -0.509$ at $(\alpha,\beta) = (0,-5)$ , clearly nonzero — the GS has converged to the wrong marginals.

Limitation. The GSC² criteria require both the joint posterior form and the full conditional densities to be available analytically. They are not applicable in fully black-box settings where Gibbs is used precisely because no analytical form exists.

CUMSUM Plot (Yu and Mykland 1994)

$\mathrm{CS}_t(X_i) = \frac{t^{-1}\sum_{j=1}^t x_i^{(j)} - \hat\mu_i}{\hat\sigma_i}, \qquad t = 1, \ldots, k$

If the chain has converged, $\mathrm{CS}_t \to 0$ smoothly. Erratic or persistent drifts indicate the chain is exploring slowly or is stuck.

Spectral Variance Estimator (Geweke 1992)

The Markov chain of sampled values is a dependent time series. The variance of the ergodic mean estimator $\hat\mu_i$ is:

$V(\hat\mu_i) = \frac{S(0)}{k} \tag{SE}$

where $S(\omega)$ is the spectral density of the series $\{x_i^{(j)}\}$ evaluated at $\omega = 0$ :

$S(\omega) = \sum_{j=-\infty}^{\infty} \gamma_j e^{-i\omega j}$

A peak in $S(\omega)$ near $\omega = 0$ indicates strong low-frequency dependence (slow mixing). In practice, $S(0)$ is estimated via a smoothed periodogram using one of three spectral windows applied to the sample autocovariances $\hat\gamma_j$ :

Window	Weight $k(j)$ for $\|j\| \leq M$
Truncated periodogram	$1$
Daniell	$\sin(\pi j/M)/(\pi j/M)$
Tukey-Hanning	$\frac{1}{2}(1 + \cos(\pi j/M))$

The bandwidth $M$ should be set to half the critical bandwidth of the spectrum — the frequency interval over which $S(\omega)$ changes appreciably. For nearly white-noise chains ( $\hat\gamma_j \approx 0$ for $j \geq 1$ ), all windows give essentially the same $\hat{S}(0)$ , and few draws are needed. For chains with strong low-frequency peaks, $M$ must be large and long chains are required.

The formula (SE) yields a confidence interval for $\hat\mu_i$ : a 95% CI is $\hat\mu_i \pm 1.96\sqrt{\hat{S}(0)/k}$ .

Why It Matters

Gibbs sampling makes Bayesian inference feasible in high-dimensional parameter spaces where direct integration is impossible.
The Chib (1995) algorithm for Marginal Data Density computation relies on Gibbs output to evaluate the posterior at a point.
Griddy-Gibbs extends the method to models with non-conjugate posteriors without requiring an importance function.
Rao-Blackwell and blocked sampling can dramatically reduce the required chain length, making large VAR estimation practical.

Open Questions

Chain mixing can fail badly for multimodal posteriors (sampler trapped near one mode); multiple restarts from dispersed starting values are a heuristic but not a guarantee.
Optimal blocking of parameters (which sets to draw jointly) has no general algorithmic solution.
For very large systems, even a single Gibbs iteration can be expensive; sequential Monte Carlo and variational approximations are increasingly used as alternatives.
SUR high-dimensional failure: In the Seemingly Unrelated Regression model with $m$ equations and $p_j$ predictors per equation, Gibbs sampling of $(\beta|\Omega, \Omega|\beta)$ requires large matrix inversions ( $\Omega^{-1} \otimes I$ ) at each step and generates highly autocorrelated draws (inefficiency factor INEF = 31.6 for off-diagonal $\Omega$ elements). Zellner and Ando (2010) show Gibbs fails to converge 92/100 trials when $p_j = 25$ , $n = 50$ . A triangular reparameterization of the SUR model yields independent normal–inverse-gamma conditionals, enabling Direct Monte Carlo (DMC) sampling with INEF = 1 and no burn-in. Percy (1992) is an early use of Gibbs specifically for the SUR predictive density problem (three-block cycle over $y_{n+1}$ , $\Phi$ , $\beta$ ), where the failure mode is less acute because $p$ (response dimension) is small. See Seemingly Unrelated Regression.