Markov Chain Monte Carlo (MCMC): Metropolis-Hastings Algorithms

Let $\{\mathbf{X}_0, \mathbf{X}_1, \mathbf{X}_2, \dots\}$ be a discrete-time, continuous-state Markov chain defined on state space $\mathcal{X} \subseteq \mathbb{R}^p$. The chain is governed by a transition kernel $P(\mathbf{x}, \mathcal{A}) = P(\mathbf{X}_{t+1} \in \mathcal{A} \mid \mathbf{X}_t = \mathbf{x})$. When the transition kernel admits a transition density $K(\mathbf{x}, \mathbf{y})$:

where $r(\mathbf{x}) = 1 - \int K(\mathbf{x}, \mathbf{y})\,d\mathbf{y}$ is the probability of remaining at $\mathbf{x}$.

1. Invariance and Detailed Balance

A probability density $\pi(\mathbf{x})$ is a stationary (invariant) distribution for the Markov chain if:

A sufficient (though not strictly necessary) condition for $\pi(\mathbf{x})$ to be stationary is the Detailed Balance (Reversibility) condition:

Proof that Detailed Balance Implies Invariance: Integrating both sides of detailed balance with respect to $\mathbf{x}$:

Since $K(\mathbf{y}, \mathbf{x})$ is a probability density integrating to 1 over $\mathcal{X}$:

Thus, detailed balance guarantees that $\pi(\mathbf{x})$ is the stationary distribution of the chain.

2. Ergodicity and the Ergodic Theorem

A Markov chain is ergodic if it is:

1. $\pi$-irreducible: Every Borel set $\mathcal{A}$ with $\pi(\mathcal{A}) > 0$ can be reached from any starting state in a finite number of steps.

2. Aperiodic: The chain does not cycle deterministically between disjoint partitions.

3. Harris recurrent: The chain revisits every non-negligible region infinitely often with probability 1.

Theorem (The Ergodic Theorem / Strong Law for MCMC): If a Markov chain $\{\mathbf{X}_t\}$ is ergodic with stationary distribution $\pi$, then for any integrable function $h: \mathcal{X} \to \mathbb{R}$ ($\mathbb{E}_\pi[|h|] < \infty$), the empirical trajectory average converges almost surely to the true expectation:

regardless of the initial starting state $\mathbf{X}_0$.

In 1953, Nicholas Metropolis, Arianna Rosenbluth, Marshall Rosenbluth, Augusta Teller, and Edward Teller published the seminal Metropolis algorithm for symmetric proposals. In 1970, W. K. Hastings generalized the algorithm to arbitrary, non-symmetric proposal distributions.

The Metropolis-Hastings Algorithm Mechanics

Let $\pi(\mathbf{x}) = \frac{\tilde{\pi}(\mathbf{x})}{Z}$ be the target probability density, where the normalizing constant $Z$ is unknown. Let $q(\mathbf{y} \mid \mathbf{x})$ be a candidate proposal distribution specifying the probability density of proposing a move from current state $\mathbf{x}$ to candidate state $\mathbf{y}$.

Given current state $\mathbf{X}_t = \mathbf{x}$:

1. Propose Candidate: Draw a candidate state $\mathbf{Y} \sim q(\mathbf{y} \mid \mathbf{x})$.

2. Compute Acceptance Probability:

Notice that the intractable normalizing constant $Z$ cancels out identically in the ratio $\frac{\tilde{\pi}(\mathbf{Y})/Z}{\tilde{\pi}(\mathbf{x})/Z}$!

3. Acceptance Step:

• Draw $U \sim \mathcal{U}(0, 1)$.
• If $U \le \alpha(\mathbf{x}, \mathbf{Y})$, accept the proposal: $\mathbf{X}_{t+1} = \mathbf{Y}$.
• Otherwise, reject the proposal: $\mathbf{X}_{t+1} = \mathbf{x}$ (the chain remains at the current state).

Proof of Detailed Balance

The effective transition density for a move from $\mathbf{x} \ne \mathbf{y}$ is:

Assume without loss of generality that $\pi(\mathbf{y}) q(\mathbf{x} \mid \mathbf{y}) \le \pi(\mathbf{x}) q(\mathbf{y} \mid \mathbf{x})$. Then:

We evaluate the forward probability flow:

We evaluate the reverse probability flow:

Both flows are strictly equal:

Detailed balance holds exactly.

The statistical performance and mixing speed of Metropolis-Hastings depends critically on the proposal density $q(\mathbf{y} \mid \mathbf{x})$.

1. Random-Walk Metropolis-Hastings

The candidate is generated by adding a zero-mean random perturbation $\mathbf{\epsilon}$ to the current state:

where $g(\mathbf{\epsilon})$ is symmetric ($g(\mathbf{\epsilon}) = g(-\mathbf{\epsilon})$), such as a Gaussian $\mathcal{N}(\mathbf{0}, \sigma_{\text{prop}}^2 \mathbf{I})$. Because $q(\mathbf{y} \mid \mathbf{x}) = g(\mathbf{y} - \mathbf{x}) = g(\mathbf{x} - \mathbf{y}) = q(\mathbf{x} \mid \mathbf{y})$, the proposal ratio cancels to 1:

This recovers the original Metropolis Algorithm.

The Proposal Scale Dilemma:

• If proposal scale $\sigma_{\text{prop}}$ is too small: Candidates $\mathbf{y}$ are close to $\mathbf{x}$, resulting in very high acceptance ($\alpha \approx 100\%$), but the chain moves via minuscule Brownian diffusion, taking millions of iterations to traverse the target distribution (high autocorrelation, poor mixing).
• If proposal scale $\sigma_{\text{prop}}$ is too large: Candidates land in low-probability tail regions where $\pi(\mathbf{y}) \ll \pi(\mathbf{x})$, resulting in near-zero acceptance ($\alpha \approx 0\%$). The chain gets stuck at the same state for thousands of iterations.

2. The Roberts-Gelman-Gilks Optimal Scaling Rule (0.234)

In 1997, Gareth Roberts, Andrew Gelman, and Walter Gilks published a foundational asymptotic theorem: Theorem (Roberts, Gelman, & Gilks, 1997): For an i.i.d. target distribution in $d$-dimensional Euclidean space $\mathbb{R}^d$ as $d \to \infty$, the optimal proposal standard deviation that maximizes the speed of diffusion of the Markov chain scales as:

Under this optimal scaling, the asymptotic average acceptance rate converges to:

For 1-dimensional targets, the optimal acceptance rate is approximately $44\%$; for $d \ge 5$ dimensions, tuning the proposal variance to maintain an acceptance rate between $20\%$ and $25\%$ provides optimal mixing.

Because MCMC samples are auto-correlated and depend on the initial arbitrary starting state $\mathbf{X}_0$, analysts must determine when the chain has reached stationarity (the burn-in period) and assess the precision of trajectory estimates.

1. Burn-In and Trace Plots

• Burn-In Period ($B$): Early iterations ($t = 1, \dots, B$) generated while the chain travels from the starting point $\mathbf{X}_0$ toward the mode of $\pi$ are discarded:

• Trace Plots: Plotting state values $X_t$ against iteration number $t$. A healthy chain exhibits a "hairy caterpillar" appearance, oscillating rapidly without secular drift.

2. Autocorrelation and Effective Sample Size (ESS)

Due to positive serial correlation $\rho_k = \mathrm{Corr}(X_t, X_{t+k})$, the variance of the sample mean is inflated by the Integrated Autocorrelation Time (IACT):

The Effective Sample Size (ESS) represents the number of equivalent independent draws:

If $N = 10,000$ and $\tau = 20$, the chain contains only $\text{ESS} = 500$ independent pieces of statistical information.

3. The Gelman-Rubin Diagnostic ($\hat{R}$)

Andrew Gelman and Donald Rubin (1992) introduced the standard multi-chain convergence metric:

• Run $M \ge 3$ parallel chains with overdispersed initial values for $N$ iterations.
• Compute the Between-Chain Variance $B / N$ and Within-Chain Variance $W$.
• The estimated marginal posterior variance is:

• The Potential Scale Reduction Factor (PSRF) is:

At convergence, between-chain variance vanishes, and $\hat{R} \to 1.0$. As a regulatory and scientific convention, convergence is diagnosed when:

In Unit 6, we demonstrated that if a Markov transition kernel $K(\mathbf{x}, \cdot)$ satisfies $\pi$-irreducibility, aperiodicity, and Harris recurrence, the distribution of state $\mathbf{X}_n$ converges in total variation distance to the stationary target distribution $\pi$:

However, Harris recurrence does not specify how rapidly this convergence occurs. If convergence is sublinear (e.g., polynomial $O(n^{-\alpha})$), an empirical MCMC run with $10^6$ iterations may remain severely contaminated by starting-value bias. To guarantee reliable simulation inference, Markov chains must satisfy Geometric Ergodicity.

1. Geometric Ergodicity and the Spectral Gap

A Markov chain with transition kernel $K$ and stationary distribution $\pi$ is geometrically ergodic if there exists a constant $\rho \in [0, 1)$ and a function $M(\mathbf{x}) < \infty$ such that:

If the bounding constant $M(\mathbf{x}) \le M$ is independent of starting state $\mathbf{x}$, the chain is uniformly ergodic, which implies exponential convergence uniform across the entire state space.

Consider the transition operator $K$ acting on the Hilbert space $L^2(\pi)$ of square-integrable functions with inner product $\langle g, h \rangle_\pi = \int g(\mathbf{x}) h(\mathbf{x}) \pi(d\mathbf{x})$. For reversible Markov chains (satisfying detailed balance), $K$ is a self-adjoint linear operator. The spectrum of $K$ is bounded within $[-1, 1]$, with the eigenvalue $1$ corresponding to constant functions.

The Absolute Spectral Gap $\gamma^*$ is defined as:

A reversible Markov chain is geometrically ergodic if and only if it possesses a positive spectral gap: $\gamma^* > 0$. The asymptotic rate of convergence is governed directly by $\rho = 1 - \gamma^*$.

2. Drift and Minorization (Foster-Lyapunov) Conditions

Proving geometric ergodicity directly from kernel definitions is notoriously difficult. Sean Meyn and Richard L. Tweedie (1993) established the Foster-Lyapunov drift condition framework, which provides practical analytical criteria for geometric ergodicity.

Theorem (Foster-Lyapunov Drift Condition): Let $C \subset \mathcal{X}$ be a small set (meaning there exists $m \ge 1$, $\epsilon > 0$, and a probability measure $\nu$ such that $K^m(\mathbf{x}, A) \ge \epsilon \nu(A)$ for all $\mathbf{x} \in C$). If there exists a Lyapunov function $V: \mathcal{X} \to [1, \infty)$ and constants $\lambda \in (0, 1)$ and $b < \infty$ such that:

then the Markov chain is geometrically ergodic.

Intuition: Outside the compact core $C$, the expected value of the Lyapunov function strictly contracts by factor $\lambda < 1$ at every step. This powerful inward drift prevents the random walk from drifting into the tail regions, guaranteeing rapid geometric convergence back to the high-density core of $\pi$.

In Random Walk Metropolis (RWM) sampling, the practitioner must select the proposal covariance matrix $\mathbf{\Sigma}_{\text{prop}} = \sigma_d^2 \mathbf{I}_d$. If the proposal step size $\sigma_d$ is chosen too small, the chain accepts almost every proposal ($\alpha \approx 100\%$), but diffuses through the state space with agonizing slowness ($O(n^{-1/2})$ random walk exploration). If $\sigma_d$ is chosen too large, the proposals land in distant tail regions with near-zero density, causing almost all steps to be rejected ($\alpha \approx 0\%$) and leaving the chain stuck at its current state.

Andrew Gelman, Gareth O. Roberts, and Walter R. Gilks (1997) proved one of the most celebrated theorems in computational statistics: finding the optimal proposal scale and acceptance rate.

1. The Asymptotic Diffusion Limit

Consider an independent and identically distributed $d$-dimensional target distribution:

Let the proposal distribution be Gaussian with scaled variance:

where $\ell > 0$ is a fixed scaling parameter, and time is accelerated by a factor of $d$: $Z_t^{(d)} = X_{\lfloor dt \rfloor, 1}$.

Theorem (Roberts, Gelman, and Gilks, 1997): As dimension $d \to \infty$, the first coordinate process $Z_t^{(d)}$ converges weakly to a continuous Langevin diffusion process $Z_t$ satisfying the Stochastic Differential Equation (SDE):

where $B_t$ is standard Brownian motion, and the diffusion speed function $h(\ell)$ is given by:

where $I = \int_{-\infty}^\infty \left( \frac{f'(x)}{f(x)} \right)^2 f(x)\,dx$ is the Fisher information of a single component.

2. Maximizing Diffusion Speed and the 0.234 Rule

Maximizing the rate of exploration of the state space corresponds to maximizing the diffusion speed coefficient $h(\ell)$ with respect to $\ell$:

Solving this optimization numerically reveals that the maximum occurs when $\ell^* = 2.38 / \sqrt{I}$. Evaluating the asymptotic acceptance rate at this optimal scaling yields:

This establishes the universal 23.4% Acceptance Rate Rule: for high-dimensional random walk Metropolis, proposal dispersion should be tuned so that approximately $23.4\%$ of proposed moves are accepted. (For 1-dimensional targets, the optimal acceptance rate is approximately $0.44$).

Because MCMC iterations are sequentially correlated, assessing convergence and computing genuine estimation precision requires specialized diagnostic metrics.

1. Effective Sample Size (ESS) and Autocorrelation

For a stationary Markov chain with variance $\sigma^2$ and lag-$k$ autocorrelation $\rho_k = \mathrm{Corr}(X_t, X_{t+k})$, the variance of the empirical sample mean $\bar{X}_n = \frac{1}{n} \sum_{t=1}^n X_t$ is:

The quantity:

is the Integrated Autocorrelation Time. It represents the number of correlated MCMC draws required to obtain the statistical information equivalent to one independent draw. The Effective Sample Size (ESS) is:

In practice, the infinite sum $\sum \rho_k$ is truncated using Geyer's initial positive sequence estimator to prevent noise accumulation from higher-order empirical autocorrelation estimates.

2. Rank-Normalized Split-$\hat{R}$

The classical Gelman-Rubin $\hat{R}$ diagnostic compares between-chain variance $B$ to within-chain variance $W$ across $M$ parallel chains. However, Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner (2021) identified two fatal flaws in classical $\hat{R}$:

1. It fails if the distribution has heavy tails or infinite variance.
2. It fails to detect within-chain non-stationarity if the chains explore different sections of the state space.

The modern Rank-Normalized Split-$\hat{R}$ solves these deficiencies:

1. Split-Chain Verification: Each chain of length $n$ is split in half into two sub-chains of length $n/2$, creating $2M$ chains. This transforms within-chain drift into between-chain variance.

2. Rank Transformation: The pooled draws from all $2M$ chains are replaced by their normalized fractional ranks:

This maps the empirical draws onto standard normal quantiles, neutralizing heavy-tail outliers.

1. Compute $\hat{R}$ on the rank-transformed data:

A value of $\hat{R} \le 1.01$ indicates satisfactory mixing and chain convergence across both tails and central mass.

All classical MCMC algorithms—including Metropolis-Hastings and Gibbs sampling—are constructed to satisfy the Detailed Balance condition:

Detailed balance guarantees reversibility: viewing the Markov trajectory forward or backward in time reveals identical statistical distributions.

However, reversibility is merely a mathematically convenient sufficient condition for stationarity, not a necessary one. The necessary and sufficient condition for $\pi$ to be invariant is the Global Balance equation:

Persi Diaconis, Susan Holmes, and Radford Neal (2000) demonstrated that enforcing reversibility actually harms simulation efficiency: reversible chains are forced to move via diffusion, repeatedly back-tracking over previously visited states. By breaking detailed balance and constructing Non-Reversible Markov Processes, one introduces directed momentum, dramatically reducing autocorrelation and accelerating convergence.

1. The Continuous-Time PDMP Framework

Modern non-reversible simulation operates via Piecewise Deterministic Markov Processes (PDMP), where state trajectories consist of continuous deterministic motion punctuated by discrete stochastic events.

Let $(\mathbf{x}(t), \mathbf{v}(t)) \in \mathbb{R}^d \times \mathbb{R}^d$ represent continuous position and discrete velocity. A PDMP is governed by three operators:

1. Deterministic Drift: Between jump events, the state evolves along deterministic vector fields: $\frac{d\mathbf{x}}{dt} = \mathbf{v}$.

2. Jump Rate (Intensity): Jump events occur as an inhomogeneous Poisson process with state-dependent event rate $\lambda(\mathbf{x}, \mathbf{v})$.

3. Transition Kernel: When a jump event occurs, the velocity updates instantaneously to $\mathbf{v}^* \sim Q((\mathbf{x}, \mathbf{v}), \cdot)$, while position remains continuous.

2. The Zig-Zag Sampler and the Bouncy Particle Sampler

Two breakthrough PDMP algorithms illustrate this non-reversible paradigm:

1. The Zig-Zag Sampler (Bierkens, Fearnhead, and Roberts 2019):

The velocity components are restricted to discrete signs $v_i \in \{-1, +1\}$. The particle travels in straight lines parallel to hypercube diagonals. Component $i$ flips velocity sign $v_i \to -v_i$ with event rate:

where $U(\mathbf{x}) = -\log \pi(\mathbf{x})$. When the particle moves uphill against the potential (away from posterior mass), the switching rate increases, flipping direction back toward the mode.

2. The Bouncy Particle Sampler (BPS) (Bouchard-Côté, Vollmer, and Doucet 2018):

Velocity moves freely on the unit sphere $\|\mathbf{v}\|_2 = 1$. When an event occurs at rate $\lambda(\mathbf{x}, \mathbf{v}) = \max(0, \langle \mathbf{v}, \nabla U(\mathbf{x}) \rangle)$, the velocity undergoes a specular reflection off the gradient hyperplane $\nabla U(\mathbf{x})$:

Because PDMPs require zero Metropolis reject-and-stay steps and accept continuous gradient information, they scale effectively to tall-data settings: Poisson thinning allows unbiased gradient evaluation using small mini-batches of data, enabling exact Bayesian MCMC on billion-row datasets.

In many scientific modeling challenges, the dimension of the parameter space itself is an unknown random variable. Examples include:

• Variable selection in regression: how many and which covariates should be included?
• Gaussian mixture modeling: how many mixture components $K \in \{1, 2, \dots, K_{\max}\}$ best explain the data?
• Change-point analysis: how many structural breaks occurred in an economic time series?

Let $\mathcal{M}_k$ denote a model indexed by $k \in \mathcal{K}$, where model $k$ possesses parameter vector $\boldsymbol{\theta}_k \in \mathbb{R}^{d_k}$. The state space of the Markov chain is the union of spaces of varying dimensions:

Standard Metropolis-Hastings fails because the target distribution has support on spaces of differing Lebesgue measures; moving between spaces of different dimensions causes the transition Jacobian to become non-square.

Peter J. Green (1995) formulated Reversible Jump Markov Chain Monte Carlo (RJMCMC), solving trans-dimensional simulation by enforcing dimension matching.

1. Dimension Matching and the Auxiliary Random Vector

Suppose the chain is at state $(k, \boldsymbol{\theta})$ of dimension $d_k$, and proposes a jump to model $(k^*, \boldsymbol{\theta}^*)$ of dimension $d_{k^*}$ (assume without loss of generality that $d_{k^*} > d_k$).

To make the jump mathematically reversible, Green introduced an auxiliary random vector $\mathbf{u} \sim q(\mathbf{u})$ of dimension $d_u$ such that the total dimensions match:

where $\mathbf{u}^*$ is an auxiliary vector generated when proposing the reverse move.

A deterministic, invertible, differentiable mapping $g$ is constructed:

Because $g$ is a diffeomorphism between spaces of equal dimension, the transformation possesses a well-defined square Jacobian matrix:

2. The Reversible Jump Acceptance Probability

Theorem (Green 1995): The detailed balance condition across models is satisfied if the acceptance probability for jumping from $(k, \boldsymbol{\theta})$ to $(k^*, \boldsymbol{\theta}^*)$ is:

where $j(k \to k^*)$ is the probability of proposing a model jump from $k$ to $k^*$.

The inclusion of the Jacobian determinant $|\det(\mathbf{J})|$ accounts precisely for the expansion or contraction of phase space volume during the dimension-changing transformation, ensuring that RJMCMC samples from the exact joint posterior distribution over both model space and parameter spaces.

When the target posterior $\pi(\mathbf{x})$ possesses multiple isolated modes separated by deep valleys of near-zero probability, standard Random Walk Metropolis and Hamiltonian Monte Carlo fail completely: the proposal probability of crossing an energy barrier decays exponentially with barrier height, trapping the chain in a single local mode for millions of iterations.

Charles J. Geyer (1991) and Koji Hukushima and Koji Nemoto (1996) introduced Parallel Tempering (also known as Replica Exchange MCMC) to traverse multimodal landscapes.

1. The Tempered Ladder System

We construct an ensemble of $K$ parallel replicas running at different "temperatures" indexed by inverse temperatures (temperatures $\beta_k \in (0, 1]$):

Each replica $k$ explores a tempered target distribution:

• At the target temperature $\beta_1 = 1$, $\pi_1(\mathbf{x}) = \pi(\mathbf{x})$ is the exact cold distribution of scientific interest.
• At high temperatures ($\beta_K \ll 1$), the potential energy barriers are flattened: high-temperature chains diffuse effortlessly between isolated modes.

2. Replica Swap Moves

In addition to standard MCMC updates within each chain, the algorithm proposes swapping states between adjacent temperature rungs $\mathbf{x}_j$ and $\mathbf{x}_{j+1}$ with Metropolis-Hastings acceptance probability:

When a high-temperature chain discovers a new mode, the swap moves propagate that state downward through the temperature ladder to the cold target chain $\beta_1 = 1$, achieving ergodic mixing across severely fractured multimodal posteriors.

A classic benchmark for diagnosing MCMC algorithms in hierarchical Bayesian modeling is Neal's Funnel:

In this 10-dimensional distribution, as $y$ becomes negative, the variance of $x_i$ shrinks exponentially, forming a narrow, high-curvature "neck" (the funnel).

When testing algorithms on $10^5$ iterations:

• Random Walk Metropolis: Becomes completely trapped in the wide opening of the funnel ($y > 0$), failing to visit the narrow neck and yielding an erroneous $\hat{R} = 1.34$.
• Standard HMC with Fixed Step Size: Experiences widespread divergent transitions as the particle bounces uncontrollably in the high-curvature neck where gradients explode.
• NUTS with Adaptive Dual Averaging and Non-Centered Parameterization: Setting $\tilde{x}_i \sim \mathcal{N}(0, 1)$ with $x_i = \tilde{x}_i e^{y/2}$ completely untangles the funnel geometry, eliminating all divergences and achieving $\hat{R} \le 1.002$ across all parameters.

Unit 6 developed the mathematical foundations of Markov Chain Monte Carlo (MCMC) algorithms. We established the conditions of $\pi$-irreducibility, aperiodicity, and Harris recurrence necessary to guarantee that the transition kernel $K(\mathbf{x}, \cdot)$ converges in total variation to the stationary distribution $\pi$. We formulated the Metropolis-Hastings algorithm, proving that enforcing the detailed balance condition $\pi(\mathbf{x}) T(\mathbf{x}, \mathbf{y}) = \pi(\mathbf{y}) T(\mathbf{y}, \mathbf{x})$ guarantees $\pi$ invariance. We analyzed Random Walk Metropolis and the Roberts-Gelman-Gilks optimal scaling theorem, which proves that in high dimensions, the optimal proposal dispersion yields an acceptance rate of exactly $23.4\%$. We analyzed geometric ergodicity through Foster-Lyapunov drift and minorization conditions, linking exponential convergence to the existence of a positive absolute spectral gap. Finally, we surveyed convergence diagnostics, detailing rank-normalized split-$\hat{R}$, effective sample size (ESS) based on integrated autocorrelation time, and parallel tempering for navigating multimodal posterior landscapes.

While the Gelman-Rubin $\hat{R}$ statistic requires running multiple independent parallel Markov chains, John Geweke (1992) proposed a spectral time-series diagnostic for evaluating the stationarity of a single MCMC chain.

Geweke's method compares the sample mean of the early portion of the chain (typically the first $10\%$, denoted $\bar{X}_A$) with the sample mean of the concluding portion of the chain (typically the last $50\%$, denoted $\bar{X}_B$), omitting the middle $40\%$ to ensure independence between the two windows. The test statistic is a two-sample asymptotic z-score:

where $\hat{S}_A(0)$ and $\hat{S}_B(0)$ are spectral density estimates of the asymptotic variance at zero frequency, computed using Parzen or Tukey-Hanning spectral windows to account for serial autocorrelation. Under the null hypothesis that the Markov chain has achieved stationarity, $Z$ is asymptotically distributed as standard normal $\mathcal{N}(0, 1)$. Values of $|Z| > 2$ provide strong statistical evidence of non-stationarity, warning the analyst that burn-in was insufficient.