Raftery-Lewis (1992) How Many Iterations in the Gibbs Sampler?

Summary

Raftery and Lewis (1992) provide a practical method for determining how many iterations a Gibbs sampler needs to estimate a given quantile $q$ of a parameter to within tolerance $r$ with probability $s$. Given a short pilot run, the method fits a two-state first-order Markov chain to the binarized output and derives the required burn-in $M$, thinning interval $k$, and total chain length $N$ analytically. It is more targeted than Gelman-Rubin: rather than assessing whether convergence has been reached, it tells you how long to run before you start.

Key Claims

• Binarize the chain at the target quantile: $Z_j = \mathbf{1}[\theta^{(j)} \leq \hat{t}]$, where $\hat{t}$ is the pilot-estimated $q$-th quantile. Under mild conditions $\{Z_j\}$ is approximately a first-order Markov chain.
• Estimate transition probabilities from pilot run: $\hat\alpha = n_{01}/(n_{00}+n_{01})$ and $\hat\beta = n_{10}/(n_{10}+n_{11})$ where $n_{ij}$ counts transitions from state $i$ to $j$.
• Test first-order adequacy with a $G^2$ statistic (deviance) comparing first-order vs. second-order models. If rejected, thin by $k$ until first-order holds.
• Required burn-in: $M \approx \lceil \log\!\bigl(\varepsilon(\hat\alpha+\hat\beta) / \max(\hat\alpha,\hat\beta)\bigr) / \log|1-\hat\alpha-\hat\beta| \rceil$ for pre-specified tolerance $\varepsilon$ on closeness to stationarity.
• Required total chain length (after burn-in, before thinning): $N_{\min} \approx \hat{N} \cdot k$ where $\hat{N} = 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 simple problems (low autocorrelation), typical recommendations are $N \approx 3{,}000$–$5{,}000$. For slowly mixing chains (high posterior correlation), $N$ can exceed $50{,}000$.
• Pilot run of $n_0 \approx 2{,}000$ draws is sufficient to estimate $\hat\alpha$, $\hat\beta$ reliably in most applications.
• The method is quantile-specific: if multiple quantiles matter, use the most demanding one (smallest $q$ or largest, closest to 0 or 1).

Concepts Introduced or Extended

• MCMC Convergence Diagnostics
• Gibbs Sampler

Entities Mentioned

• None requiring wiki pages

Quotes

"The method described here is based on the situation where there is a single long run of the Gibbs sampler, as practiced by Geman and Geman (1984) and Besag, York and Mollié (1991), for example."

"Sometimes, although not always, this may be exploited in the first way by reducing the value of $q$, to obtain the same result with less total iteration."

My Take

The method's main strength is that it gives a calibrated answer tied to a specific inferential goal: if you care about the 2.5th percentile to $\pm$0.005, the output tells you exactly how long to run and how long to burn in. The first-order Markov assumption on the binarized chain is testable (via the $G^2$ comparison) and thinning rescues it when violated. The main limitation is that the method is univariate: it applies one parameter at a time and ignores joint behavior. In high-dimensional models with correlated parameters, the worst-case parameter determines the chain length. In practice, the method is best used as a lower bound — run at least this many iterations — and combined with Gelman-Rubin multi-chain diagnostics for a fuller assessment.