1. Purpose of Iterative Simulation
In Bayesian analysis, our goal is to summarize the posterior distribution
$p(θ \mid y)$ of the model parameters.
When direct analytical integration is impossible, we use Markov chain Monte Carlo (MCMC)—for example, the Metropolis or Gibbs algorithms—to generate a sequence of dependent random draws:
$θ^{(1)}, θ^{(2)}, \dots, θ^{(S)} \;\sim\; p(θ \mid y)$
From these draws, we compute:
- Posterior summaries: mean, median, variance, quantiles.
- Posterior intervals: e.g., 50% or 95% intervals from percentiles.
- Functions of parameters: any scalar estimand ψ = h(θ).
- Posterior predictive distributions: simulate future or missing data,
$\tilde{y}^{(s)} \sim p(\tilde{y} \mid θ^{(s)})$.
So, the simulated sample from $p(θ\mid y)$ becomes the foundation for all inference.
2. Challenges Unique to Iterative (MCMC) Simulation
Compared to independent Monte Carlo draws, MCMC brings two extra difficulties.
(a) Insufficient iterations → non-convergence
If we stop too early, the chain has not yet reached the stationary distribution.
Early iterations are still influenced by the starting values,
so the draws may be biased toward initial conditions.
Above images are illustrate this:
early samples cluster near starting points, not across the true posterior region.
(b) Autocorrelation between successive draws
Even after convergence, successive samples are correlated, not independent.
This correlation inflates Monte Carlo error and reduces efficiency.
Visually, a scatterplot of consecutive draws shows “patchy” or “banded” clustering—evidence of slow movement of the chain.
Thus, the effective sample size (ESS) is smaller than the total S.
An approximate ESS formula appears later.
3. General Strategy to Handle These Problems
The propose three integrated strategies:
- Run multiple chains with overdispersed starting points
→ to evaluate whether chains converge to the same region. - Monitor convergence quantitatively, comparing variation
within each chain vs. between chains.
→ When these two sources of variation become similar, convergence is likely. - Improve efficiency if needed by modifying the algorithm:
- Reparameterize the model
- Tune step sizes, proposal scales, or update blocks
- Or use more advanced methods (Hamiltonian Monte Carlo, etc.)
4. Warm-up (formerly “burn-in”)
To reduce the effect of the starting points, we discard the early part of each simulated chain.
- The authors call this warm-up, avoiding the older term “burn-in.”
- Default rule: discard roughly the first half of each sequence.
- Example: run 200 iterations → discard the first 100 and keep the last 100.
- If convergence not reached, extend the run and again discard all early draws.
- The fraction discarded can vary by problem; for example,
a simple Gibbs sampler might need only a few iterations.
5. Thinning (skipping iterations)
After convergence, one may consider keeping every k-th draw (e.g., 1 of 10) to reduce autocorrelation.
- Thinning is not necessary for valid inference; it only helps if
storage or memory limits are severe. - The authors recommend setting k just large enough to limit total stored draws to about 1000 if parameters are numerous.
- Otherwise, keep all post-warm-up draws.
6. Multiple Sequences with Overdispersed Starting Points
To assess convergence properly, it is essential to run multiple independent chains.
- Start each chain from widely separated (overdispersed) points in parameter space.
- These starting values can come from crude estimates
or from more sophisticated approximations.
- These starting values can come from crude estimates
Reason:
If chains starting far apart all end up exploring the same region of the posterior, we can be confident that the Markov chains have mixed well and reached stationarity.
Before convergence: large between-chain differences.
After convergence: chains overlap; between- and within-variances equalize.
7. Monitoring Scalar Estimands (ψ)
Convergence must be checked separately for each scalar quantity of interest ψ.
These can include:
- Individual parameters $θ_j$
- Derived quantities (ratios, differences, predictions)
- Log-posterior density $\log p(θ\mid y)$
Monitoring each ψ separately helps detect local convergence issues even if the overall likelihood looks stable.
8. Understanding “Mixing” and “Stationarity”
Above image is highlights two distinct convergence requirements:
- Mixing:
Chains must explore the same region of the parameter space.- Two separate chains can each look stable individually yet be centered at different means.
- Only by comparing them together can we see they haven’t mixed.
- Stationarity:
Each individual chain must have stabilized over time.- A chain might wander systematically, even if multiple chains overlap.
Hence, both between-sequence mixing and within-sequence stationarity must be checked.
9. Splitting Chains to Diagnose Convergence
To assess both criteria simultaneously, we use the split-chain method (the basis for modern $\hat{R}$).
Steps:
- Start with m independent chains (after warm-up).
Example: 5 chains, each 1000 iterations. - Discard first ½ of each for warm-up → remaining 500 each.
- Split each remaining chain into two halves → 10 chains of 250 each.
- $m = 10, n = 250$
- Treat these as independent sequences for the calculations below.
10. Between- and Within-Sequence Variances
For each scalar estimand ψ:
Let $ψ_{ij}$ be the i-th draw of the j-th chain (after splitting).
$\begin{aligned} ψ_{\cdot j} &= \frac{1}{n}\sum_{i=1}^{n} ψ_{ij} \quad &\text{(mean of chain j)}\\ ψ_{\cdot\cdot} &= \frac{1}{m}\sum_{j=1}^{m} ψ_{\cdot j} \quad &\text{(grand mean)}\\[4pt] W &= \frac{1}{m}\sum_{j=1}^{m} s_j^2, \quad s_j^2 = \frac{1}{n-1}\sum_{i=1}^{n} (ψ_{ij} – ψ_{\cdot j})^2 \quad &\text{(within-chain variance)}\\[4pt] B &= \frac{n}{m-1}\sum_{j=1}^{m} (ψ_{\cdot j} – ψ_{\cdot\cdot})^2 \quad &\text{(between-chain variance)} \end{aligned}$
11. Estimating the Posterior Variance
The overall posterior variance of ψ is estimated by combining B and W:
$\mathrm{Var}^+(ψ\mid y) = \frac{n-1}{n}W + \frac{1}{n}B$
- $W$ underestimates the true variance because individual chains haven’t fully explored the distribution.
- $B$ overestimates variance if starting points were too dispersed.
- The weighted average $\mathrm{Var}^+$ is unbiased under stationarity and conservative otherwise.
12. Potential Scale Reduction Factor ($\hat R$)
To quantify how much the scale of the current simulated distribution could shrink with infinite iterations, compute:
$\hat R = \sqrt{\frac{\mathrm{Var}^+(ψ\mid y)}{W}}$
Interpretation
- $\hat R \to 1$ as $n \to \infty$ → convergence.
- $\hat R > 1$ → chains still differ; more simulation needed.
- In practice:
- $\hat R < 1.05$ = acceptable;
- $\hat R \approx 1.00\!-\!1.01$ = excellent.
The diagnostic measures how much better we could estimate the variance if chains ran forever—hence “potential scale reduction.”
13. Worked Example – Bivariate Normal via Metropolis
The authors demonstrate the method with the bivariate standard normal target:
- Target: $p(θ \mid y) = N(0, I_2)$
- Sampler: Metropolis algorithm with small step size (deliberately inefficient).
- Parallel runs: 5 independent sequences, each run for 50, 500, 2000, 5000 iterations.
Results
→ After about 2000–5000 iterations, the chains converge adequately;
$\hat R$ values approach 1, and posterior intervals match the true normal distribution.
Even with perfect convergence, posterior estimates still vary slightly because Monte Carlo noise remains; thousands of draws are often required for precise tail estimates.
14. Advantages of the $\hat R$ Diagnostic
- No need for time-series plots.
Visual inspection is subjective and unreliable, especially with many parameters. - Scalable.
Works automatically for large hierarchical models with hundreds of monitored quantities. - Statistical basis.
Built on variance decomposition (within + between).
However, since it relies on means and variances, the method works best when marginal posteriors are approximately normal.
For parameters with high skewness, multimodality, or extreme quantiles, also inspect chain behavior for those tail areas separately.
15. Summary of Best Practices
| Step | Purpose | Typical Action |
|---|---|---|
| 1. Multiple chains | Assess global convergence | ≥ 2 chains, widely separated starts |
| 2. Warm-up | Remove early bias | Discard ~ ½ of each chain |
| 3. Split chains | Check stationarity + mixing | Split remaining chains in half |
| 4. Compute W, B, $\hat R$ | Quantitative convergence test | Aim for $\hat R ≈ 1.00–1.05$ |
| 5. ESS | Gauge simulation precision | Use later formula; higher = better |
| 6. Use all post-warm-up draws | Posterior inference | Don’t thin unless storage is limiting |
| 7. Report diagnostics | Transparency | Include $\hat R$ and ESS in results |
Final summary:
In Bayesian inference via MCMC, reliable conclusions depend on verifying that simulated chains have converged.
This is achieved by running multiple overdispersed chains, discarding early warm-up iterations, splitting chains to test both mixing and stationarity, and computing the potential scale reduction factor $\hat R$.
When $\hat R ≈ 1$, the draws can be treated as samples from $p(θ | y)$, allowing accurate posterior summaries and predictive inference.