The Gibbs Sampler & Bayesian Computational Inference

Let $\mathbf{\theta} = (\theta_1, \theta_2, \dots, \theta_p)^T \in \Theta$ be a parameter vector with joint posterior distribution $\pi(\mathbf{\theta} \mid \mathbf{y})$. Let $\mathbf{\theta}_{-j} = (\theta_1, \dots, \theta_{j-1}, \theta_{j+1}, \dots, \theta_p)^T$ denote the subvector of all parameters excluding $\theta_j$.

The Full Conditional Distribution of parameter $\theta_j$ is the conditional distribution of $\theta_j$ given all remaining parameters and the data:

1. The Systematic Scan Gibbs Sampler

Given an initial state $\mathbf{\theta}^{(0)} = (\theta_1^{(0)}, \dots, \theta_p^{(0)})^T$, iteration $t+1$ executes a deterministic sweep through all $p$ coordinates:

1. Sample $\theta_1^{(t+1)} \sim \pi(\theta_1 \mid \theta_2^{(t)}, \theta_3^{(t)}, \dots, \theta_p^{(t)}, \mathbf{y})$.
2. Sample $\theta_2^{(t+1)} \sim \pi(\theta_2 \mid \theta_1^{(t+1)}, \theta_3^{(t)}, \dots, \theta_p^{(t)}, \mathbf{y})$.

$\dots$ $j$. Sample $\theta_j^{(t+1)} \sim \pi(\theta_j \mid \theta_1^{(t+1)}, \dots, \theta_{j-1}^{(t+1)}, \theta_{j+1}^{(t)}, \dots, \theta_p^{(t)}, \mathbf{y})$. $\dots$ $p$. Sample $\theta_p^{(t+1)} \sim \pi(\theta_p \mid \theta_1^{(t+1)}, \dots, \theta_{p-1}^{(t+1)}, \mathbf{y})$.

Notice that as each new parameter $\theta_j^{(t+1)}$ is drawn, it is immediately used in the conditioning sets of all subsequent updates within the same iteration.

2. The Gibbs Sampler as Metropolis-Hastings with Acceptance $\alpha \equiv 1$

A profound theoretical insight is that the Gibbs sampler is mathematically equivalent to a multi-stage Metropolis-Hastings algorithm where the acceptance probability is identically 100%:

Theorem: Consider updating coordinate $j$ using candidate proposal $q(\mathbf{\theta}^* \mid \mathbf{\theta}) = \pi(\theta_j^* \mid \mathbf{\theta}_{-j}, \mathbf{y}) \cdot I(\mathbf{\theta}_{-j}^* = \mathbf{\theta}_{-j})$. The Metropolis-Hastings acceptance ratio is:

By definition of conditional probability:

Substituting into the acceptance ratio:

Every single candidate draw is accepted with probability $\alpha \equiv 1.0$! Zero proposals are rejected.

In many statistical models, the observed data likelihood $p(\mathbf{y} \mid \mathbf{\theta})$ is analytically intractable, but the complete data likelihood $p(\mathbf{y}, \mathbf{z} \mid \mathbf{\theta})$ involving unobserved latent variables $\mathbf{z}$ is simple and conjugate. In 1987, Martin Tanner and Wing Hung Wong formulated the Data Augmentation Framework.

The Tanner-Wong Two-Step Iteration

1. Imputation Step (I-Step): Sample the missing latent variables from their predictive distribution given current parameters:

2. Posterior Step (P-Step): Sample the parameters from their complete-data posterior distribution given augmented data:

As $t \to \infty$, the marginal sequence $\{\mathbf{\theta}^{(t)}\}$ converges in distribution to the exact observed posterior $\pi(\mathbf{\theta} \mid \mathbf{y})$.

The Albert-Chib Algorithm for Bayesian Probit Regression

Consider a binary probit regression model where observed responses $y_i \in \{0, 1\}$ follow:

The likelihood function $\prod_{i=1}^n \Phi(\mathbf{x}_i^T \mathbf{\beta})^{y_i} [1 - \Phi(\mathbf{x}_i^T \mathbf{\beta})]^{1 - y_i}$ is non-conjugate with respect to a Gaussian prior $\mathbf{\beta} \sim \mathcal{N}(\mathbf{\beta}_0, \mathbf{\Sigma}_0)$.

In 1993, James Albert and Siddhartha Chib introduced continuous latent variables $Z_i \sim \mathcal{N}(\mathbf{x}_i^T \mathbf{\beta}, 1)$ such that:

Under this formulation, the full conditionals are analytically tractable standard distributions:

1. Update Latent Variables: For each $i = 1, \dots, n$:

• If $y_i = 1$, sample $Z_i \sim \mathcal{N}(\mathbf{x}_i^T \mathbf{\beta}, 1)$ truncated to the positive interval $(0, \infty)$.
• If $y_i = 0$, sample $Z_i \sim \mathcal{N}(\mathbf{x}_i^T \mathbf{\beta}, 1)$ truncated to the negative interval $(-\infty, 0]$.

2. Update Regression Coefficients: Conditional on continuous vector $\mathbf{Z} = (Z_1, \dots, Z_n)^T$, the model is standard linear regression $\mathbf{Z} = \mathbf{X}\mathbf{\beta} + \mathbf{\epsilon}$ with $\mathbf{\epsilon} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. With prior $\mathbf{\beta} \sim \mathcal{N}(\mathbf{\beta}_0, \mathbf{\Sigma}_0)$:

The Albert-Chib sampler transforms non-linear categorical inference into linear Gaussian updates.

Consider the classical multiple linear regression model:

where $\mathbf{y} \in \mathbb{R}^n$ is the response vector, $\mathbf{X} \in \mathbb{R}^{n \times p}$ is the design matrix, $\mathbf{\beta} \in \mathbb{R}^p$ is the regression coefficient vector, and $\sigma^2 > 0$ is the residual variance.

We place independent semi-conjugate priors on $\mathbf{\beta}$ and $\sigma^2$:

• $\mathbf{\beta} \sim \mathcal{N}_p(\mathbf{\beta}_0, \mathbf{\Sigma}_0)$ (Gaussian prior).
• $\sigma^2 \sim \text{Inverse-Gamma}(a_0, b_0)$ with density $p(\sigma^2) \propto (\sigma^2)^{-(a_0 + 1)} \exp(-b_0 / \sigma^2)$.

The Full Conditional Posterior Distributions

1. Full Conditional of $\mathbf{\beta}$ given $\sigma^2$ and $\mathbf{y}$:

where:

2. Full Conditional of $\sigma^2$ given $\mathbf{\beta}$ and $\mathbf{y}$:

where:

The Gibbs sampler alternates between drawing $\mathbf{\beta}^{(t+1)} \sim \mathcal{N}_p(\mathbf{\beta}_n, \mathbf{\Sigma}_n)$ using Cholesky factorization and drawing $\sigma^{2(t+1)} \sim \text{Inverse-Gamma}(a_n, b_n)$ via inverted Gamma sampling.

In high-dimensional parameter spaces ($d \ge 50$), Random Walk Metropolis and Gibbs sampling suffer from diffusive behavior: the distance traveled scales only as $\sqrt{n}$. To escape diffusion, Radford Neal (2011) adapted physical mechanics into Hamiltonian Monte Carlo (HMC), originally introduced by Simon Duane, Anthony Kennedy, Brian Pendleton, and Duncan Roweth in 1987.

HMC translates Bayesian posterior sampling into the physical simulation of a particle sliding frictionless over a multi-dimensional energy potential.

1. The Hamiltonian Dynamical System

Let $\mathbf{q} \in \mathbb{R}^d$ represent the position vector (the parameters of interest), with potential energy equal to the negative unnormalized log-posterior:

We introduce fictitious auxiliary momentum variables $\mathbf{p} \in \mathbb{R}^d$, typically drawn from a multivariate Gaussian distribution $\mathbf{p} \sim \mathcal{N}(\mathbf{0}, \mathbf{M})$, where $\mathbf{M}$ is a symmetric positive-definite mass matrix. The kinetic energy is:

The total energy of the physical system is the Hamiltonian:

The time evolution of the system follows Hamilton's equations of motion:

Three fundamental physical properties guarantee the validity of HMC:

1. Conservation of Energy: $\frac{dH}{dt} = \sum_i \left( \frac{\partial H}{\partial q_i} \frac{dq_i}{dt} + \frac{\partial H}{\partial p_i} \frac{dp_i}{dt} \right) = \sum_i \left( \frac{\partial H}{\partial q_i} \frac{\partial H}{\partial p_i} - \frac{\partial H}{\partial p_i} \frac{\partial H}{\partial q_i} \right) = 0$.

2. Reversibility: Negating the momentum $\mathbf{p} \to -\mathbf{p}$ reverses the trajectory precisely.

3. Volume Preservation (Liouville's Theorem): The phase space flow $(\mathbf{q}(0), \mathbf{p}(0)) \mapsto (\mathbf{q}(t), \mathbf{p}(t))$ has a Jacobian determinant identically equal to $1.0$, meaning phase space volume is strictly conserved: $|\det(J)| = 1$.

2. The Symplectic Leapfrog Integrator

Because continuous differential equations cannot be solved analytically for arbitrary posteriors, HMC discretizes time into steps of size $\epsilon$ using the Leapfrog Integrator:

Unlike naive Euler integration (which spirals away due to numerical energy drift), the leapfrog integrator is symplectic: it exactly preserves an approximate "shadow Hamiltonian" $\tilde{H} = H + O(\epsilon^2)$.

After $L$ leapfrog steps from $(\mathbf{q}, \mathbf{p})$ to proposed state $(\mathbf{q}^*, \mathbf{p}^*)$, the Metropolis acceptance probability is:

Because $\Delta H = O(\epsilon^2)$ is purely due to discretization error, $\alpha$ remains close to $1.0$ (typically tuned to $65\% - 80\%$) even across vast distances in hundreds of dimensions.

While HMC overcomes random walk diffusion, its performance is sensitive to two hyperparameters: the step size $\epsilon$ and the trajectory length $L$. If $L$ is too small, the algorithm degenerates into random walk exploration. If $L$ is too large, the simulated particle loops back along its orbital trajectory, wasting expensive gradient evaluations.

Matthew D. Hoffman and Andrew Gelman (2014) introduced the No-U-Turn Sampler (NUTS) to automate trajectory termination.

1. The U-Turn Criterion

The trajectory should be terminated the instant the particle begins folding back on itself, meaning the Euclidean distance between the current position $\mathbf{q}$ and the initial position $\mathbf{q}_0$ stops increasing:

2. Recursive Binary Tree Building

To preserve exact time-reversibility and detailed balance, NUTS builds a balanced binary tree by simulating leapfrog steps forward and backward in time with equal probability. If the sub-trajectories in either direction satisfy the U-turn condition:

where $\mathbf{q}^+$ and $\mathbf{q}^-$ are the leftmost and rightmost endpoints of the tree, trajectory expansion halts. A point is then sampled from the tree trajectory using progressive multinomial sampling.

3. Dual Averaging Step Size Adaptation

NUTS tunes step size $\epsilon$ during warmup using Yurii Nesterov's primal-dual subgradient method (Dual Averaging). The step size is adjusted at each iteration to match a target acceptance rate $\delta_{\text{target}}$ (typically $0.80$ in Stan and PyMC):

This automated adaptation eliminates manual tuning, making NUTS the core computational engine of modern probabilistic programming.

Introduced by Radford M. Neal (2003), Slice Sampling provides an auxiliary variable MCMC method that samples from continuous distributions without requiring explicit proposal distribution tuning.

1. Univariate Slice Sampling

Let $f(x) \propto h(x)$ be an unnormalized target density. Introduce an auxiliary vertical variable $y \in \mathbb{R}$:

The marginal distribution of $x$ is:

Univariate slice sampling alternates between two steps:

1. Vertical Slice Draw: Given current state $x_0$, draw $y \sim \mathcal{U}(0, h(x_0))$. This defines a horizontal slice through the density:

2. Horizontal Coordinate Draw: Sample a new position $x_1$ uniformly from the horizontal slice set $S$.

2. Stepping-Out and Shrinkage Procedures

Because the exact boundaries of the slice $S$ cannot be determined analytically for general densities, Neal developed the stepping-out and shrinkage algorithms:

1. Stepping-Out: An initial interval $[L, R]$ of width $w$ is positioned randomly around $x_0$: $L = x_0 - U \cdot w$ and $R = L + w$ with $U \sim \mathcal{U}(0, 1)$. The interval is expanded outward in steps of size $w$ until both ends fall outside the slice ($h(L) < y$ and $h(R) < y$).

2. Shrinkage: A candidate point $x^*$ is proposed uniformly from $[L, R]$. If $h(x^*) \ge y$, $x^*$ is accepted. If $h(x^*) < y$, the candidate falls outside the slice; the interval $[L, R]$ is shrunk directly to $x^*$ (setting $L = x^*$ if $x^* < x_0$, or $R = x^*$ if $x^* > x_0$), and a new proposal is drawn from the narrowed interval.

Because rejected points immediately shrink the candidate interval, slice sampling adapts automatically to local curvature: the step size expands rapidly in flat plateaus and contracts automatically in steep ravines, achieving near-$100\%$ efficiency without tuning parameters.

In modern computational statistics and machine learning, Bayesian posterior approximation is divided into two competing paradigms: **Markov Chain Monte Carlo (MCMC)** and **Variational Inference (VI)**. While MCMC approaches posterior estimation via stochastic simulation, VI recasts posterior inference as a deterministic optimization problem. ### 1. The Evidence Lower Bound (ELBO) Let $\mathbf{x}$ be observed data and $\mathbf{z}$ be latent variables or model parameters with joint density $p(\mathbf{x}, \mathbf{z}) = p(\mathbf{x} \mid \mathbf{z}) p(\mathbf{z})$. We specify a family of tractable approximate distributions $\mathcal{Q} = \{ q_\phi(\mathbf{z}) \}$ and search for the variational parameters $\phi^*$ that minimize the Kullback-Leibler (KL) divergence to the exact posterior: $$\phi^* = \arg\min_\phi \mathcal{D}_{\text{KL}}(q_\phi(\mathbf{z}) \parallel p(\mathbf{z} \mid \mathbf{x})) = \arg\min_\phi \int q_\phi(\mathbf{z}) \log \left( \frac{q_\phi(\mathbf{z})}{p(\mathbf{z} \mid \mathbf{x})} \right) d\mathbf{z}$$ Because the posterior $p(\mathbf{z} \mid \mathbf{x}) = p(\mathbf{x}, \mathbf{z}) / p(\mathbf{x})$ contains the intractable evidence $p(\mathbf{x})$, the KL divergence cannot be minimized directly. However, decomposing the marginal log-likelihood reveals: $$\log p(\mathbf{x}) = \int q_\phi(\mathbf{z}) \log p(\mathbf{x})\,d\mathbf{z} = \int q_\phi(\mathbf{z}) \log \left( \frac{p(\mathbf{x}, \mathbf{z})}{q_\phi(\mathbf{z})} \cdot \frac{q_\phi(\mathbf{z})}{p(\mathbf{z} \mid \mathbf{x})} \right) d\mathbf{z}$$ $$\log p(\mathbf{x}) = \underbrace{\mathbb{E}_{q_\phi}\left[ \log \frac{p(\mathbf{x}, \mathbf{z})}{q_\phi(\mathbf{z})} \right]}_{\text{ELBO}(\phi)} + \underbrace{\mathcal{D}_{\text{KL}}(q_\phi(\mathbf{z}) \parallel p(\mathbf{z} \mid \mathbf{x}))}_{\ge 0}$$ Because the KL divergence is strictly non-negative, the **Evidence Lower Bound (ELBO)** forms a rigorous lower bound on the marginal likelihood: $$\log p(\mathbf{x}) \ge \text{ELBO}(\phi) = \mathbb{E}_{q_\phi}[\log p(\mathbf{x}, \mathbf{z})] - \mathbb{E}_{q_\phi}[\log q_\phi(\mathbf{z})]$$ Maximizing $\text{ELBO}(\phi)$ with respect to $\phi$ is mathematically identical to minimizing $\mathcal{D}_{\text{KL}}(q_\phi \parallel p(\cdot \mid \mathbf{x}))$. ### 2. Methodological Trade-Off Comparison In modern workflows, VI is routinely deployed for rapid exploratory data analysis, initial hyperparameter screening, and initializing MCMC warmup trajectories, while Hamiltonian Monte Carlo is employed for final regulatory and scientific hypothesis reporting.

In standard MCMC simulation, the practitioner faces an insurmountable theoretical obstacle: the Markov chain only converges to the stationary distribution $\pi$ asymptotically as $n \to \infty$. In any finite simulation of $N$ steps, the distribution of state $X_N$ is only approximately equal to $\pi$, and determining the exact warmup time required to render the remaining bias negligible is notoriously difficult.

James G. Propp and David B. Wilson (1996) shocked the probability community by inventing Coupling from the Past (CFTP), an algorithm that generates exact, unbiased, perfect draws from the stationary distribution $\pi$ of a Markov chain in finite computational time.

1. Monotonicity and the Coalescence Criterion

Consider a discrete or continuous state space $\mathcal{X}$ equipped with a partial order $\preceq$, with unique minimum element $\hat{0}$ and maximum element $\hat{1}$ (such as the Ising model of ferromagnetic spins on a lattice).

Suppose the Markov transition update can be represented by a deterministic update function $\phi: \mathcal{X} \times \mathcal{U} \to \mathcal{X}$ driven by independent uniform random numbers $U_t \sim \mathcal{U}$:

The update function is said to be monotone if:

2. Why Running Forward Fails vs. Coupling From the Past

If one starts two Markov chains at time $t = 0$ from different states and runs them forward until they meet (coalesce) at time $T$, the distribution at coalescence time $T$ is biased: stopping at the random time $T$ selects for trajectories that coalesce rapidly, violating stationarity.

Propp and Wilson realized that the coalescence must occur at fixed time $t = 0$ by starting chains in the infinite past:

1. Fix observation time $t = 0$.
2. Choose a starting time in the past $T_1 = -1$.
3. Run chains starting from all possible states $x \in \mathcal{X}$ from time $T_1$ up to $t = 0$, using the random sequence $(U_{-1})$.
4. If all trajectories coalesce into a single state at $t = 0$, output that common state $X_0$.
5. If trajectories have not coalesced, double the starting horizon: $T_2 = -2, -4, -8, \dots, -2^k$.
6. Crucial Coupling Principle: When restarting from $T_{k+1} = -2^{k+1}$, one must reuse the exact same random numbers $U_t$ that were used previously for $t \in [-2^k, 0)$. Only new random numbers are drawn for the extended interval $[-2^{k+1}, -2^k)$.

The Sandwiching Property: Thanks to monotonicity, one does not need to simulate trajectories from all $2^N$ possible states: one only simulates the two extreme chains starting from $\hat{0}$ and $\hat{1}$. If the top chain and bottom chain coalesce by time $t = 0$:

then every intermediate chain is squeezed between them and must also have coalesced to $X_0$.

Theorem (Propp-Wilson): The coalesced state $X_0$ is distributed exactly according to the stationary measure $\pi$, with zero burn-in bias and zero approximation error.

When fitting Bayesian models to massive datasets containing $N = 10^7$ or $10^8$ observations, standard Hamiltonian Monte Carlo and Metropolis algorithms become computationally prohibitive: computing the gradient of the log-posterior requires evaluating $N$ terms at every leapfrog step:

Max Welling and Yee Whye Teh (2011) formulated Stochastic Gradient Langevin Dynamics (SGLD), integrating stochastic gradient descent with continuous Langevin diffusion to enable scalable Bayesian MCMC on big data.

1. Mini-Batch Gradient Approximations

At step $t$, SGLD computes an unbiased stochastic estimate of the log-likelihood gradient using a small random mini-batch $S_t \subset \{1, \dots, N\}$ of size $n \ll N$ (e.g., $n = 256$):

The parameter vector updates according to:

where $\epsilon_t$ is a decaying learning rate schedule satisfying the Robbins-Monro conditions:

2. Transition from Optimization to Posterior Sampling

In the early iterations where $\epsilon_t$ is large, the gradient drift term $\frac{\epsilon_t}{2} \mathbf{g}_t$ dominates the injected Gaussian noise $\sqrt{\epsilon_t}$, causing the algorithm to act like mini-batch SGD and rapidly descend toward high-density posterior modes.

As $t \to \infty$ and $\epsilon_t \to 0$, the mini-batch sampling noise (which scales as $\epsilon_t^2$) decays much faster than the injected thermal noise $\sqrt{\epsilon_t}$. The injected noise dominates the sub-sampling error, causing the discrete trajectory to converge weakly to the continuous Langevin diffusion whose stationary distribution is the true Bayesian posterior. SGLD eliminates the costly $O(N)$ Metropolis accept-reject step entirely, enabling Bayesian inference on big data streams.

Latent Dirichlet Allocation (LDA) models the semantic composition of text corpora by assigning words to latent topics. Fitting an LDA model with $K = 50$ topics to a corpus of 100,000 documents:

• Collapsed Gibbs Sampling (Griffiths & Steyvers): Iterates sequentially through every word token in the corpus. While achieving high posterior perplexity and capturing multi-word semantic subtleties, each sweep requires 45 minutes of computation, taking 15 hours to converge over 20 sweeps.
• Stochastic Variational Inference (Hoffman et al.): Optimizes the variational ELBO using mini-batches of 256 documents. SVI reaches 95% of the optimal perplexity in only 12 minutes—a 75-fold reduction in wall-clock time.

This empirical comparison highlights the complementary roles of VI (for rapid discovery and real-time processing) and Gibbs/HMC (for high-precision scientific parameter quantification).

Unit 7 focused on advanced MCMC architectures and modern Hamiltonian dynamics. We analyzed the Gibbs sampler as an exact Metropolis-Hastings scheme with acceptance probability identically equal to 1.0, examining systematic-scan and random-scan variations. We formulated Hamiltonian Monte Carlo (HMC), demonstrating how auxiliary Gaussian momentum variables transform posterior sampling into frictionless physical motion governed by Hamilton's equations. We analyzed the volume-preserving, time-reversible properties of the symplectic Leapfrog integrator and the No-U-Turn Sampler (NUTS), which automates trajectory termination via recursive binary trees and dual averaging step size adaptation. We studied Radford Neal's Slice Sampling with stepping-out and shrinkage procedures, which adaptively navigates varying posterior scales without parameter tuning. Finally, we contrasted MCMC against Variational Inference (VI) and Stochastic Gradient Langevin Dynamics (SGLD), and surveyed Propp and Wilson's Coupling from the Past (CFTP) for exact, unbiased perfect simulation.

When developing custom MCMC samplers or diagnosing convergence failures in probabilistic programming models, following a systematic troubleshooting hierarchy prevents widespread inference errors:

1. Analytical Test Verification via Prior Simulation: Before fitting empirical data, simulate synthetic datasets from known parameter priors $\boldsymbol{\theta}_0 \sim p(\boldsymbol{\theta})$ and $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta}_0)$. Verify that the posterior credibility intervals cover the true parameters $\boldsymbol{\theta}_0$ at nominal theoretical frequencies across 1,000 independent simulation replications (Simulation-Based Calibration).

2. Detection of Geometric Funnels via Divergent Transitions: In Hamiltonian Monte Carlo, any single divergent transition signifies that the simulated trajectory has encountered a region of infinite curvature where numerical leapfrog integration breaks down. Reparameterizing the model from a centered formulation $\theta \sim \mathcal{N}(\mu, \sigma^2)$ to a non-centered formulation $\tilde{\theta} \sim \mathcal{N}(0, 1)$ with $\theta = \mu + \sigma \tilde{\theta}$ straightens the target geometry and eliminates divergences.

3. Energy Bayesian Fraction of Missing Information (E-BFMI): In HMC, compute the energy metric:

A value of $\text{E-BFMI} < 0.3$ indicates that the momentum distribution cannot effectively explore the energy state space, diagnosing poor momentum mass matrix adaptation or heavy tail pathologies.

In structural bioinformatics, portfolio optimization with cardinality constraints, and neural architecture search, objective loss functions $f(\mathbf{x})$ exhibit complex non-convex landscapes with thousands of deceptive local minima. Deterministic gradient descent gets trapped permanently in the nearest local basin of attraction.

Scott Kirkpatrick, C. Daniel Gelatt, and Mario P. Vecchi (1983) adapted the Metropolis MCMC algorithm into Simulated Annealing (SA), treating optimization as the physical process of slowly cooling molten metal to find its crystalline ground state.

At cooling iteration $k$, a proposed move $\mathbf{x}_{\text{prop}} \sim q(\cdot \mid \mathbf{x}_k)$ is accepted with probability:

When the temperature $T_k$ is high, the algorithm readily accepts uphill moves ($f(\mathbf{x}_{\text{prop}}) > f(\mathbf{x}_k)$), freely exploring the global state space. As temperature $T_k \to 0$ according to a logarithmic cooling schedule $T_k = T_0 / \log(k + 1)$ (Hajek 1988), the chain freezes into the deepest global minimum with probability approaching 1.0.