7.2 Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm is a powerful sampling technique in Bayesian statistics. It allows us to explore complex posterior distributions when direct sampling is impossible. This method is crucial for estimating parameters and making inferences in models with intractable posteriors.

The algorithm works by generating a Markov chain that converges to the target distribution. It uses a proposal distribution to suggest new states and an acceptance probability to decide whether to move. This process enables efficient exploration of parameter spaces in various Bayesian models.

Overview of Metropolis-Hastings algorithm

• Fundamental sampling technique in Bayesian statistics enables exploration of complex posterior distributions
• Iterative algorithm generates samples from target distribution without direct sampling
• Crucial for estimating parameters and making inferences in Bayesian models with intractable posteriors

Markov Chain Monte Carlo

MCMC fundamentals

• Sampling method constructs Markov chain with desired stationary distribution
• Relies on ergodic property ensures chain converges to target distribution regardless of starting point
• Generates dependent samples through iterative process
• Utilizes random walks to explore parameter space efficiently

Role in Bayesian inference

• Facilitates sampling from posterior distributions in complex Bayesian models
• Enables estimation of marginal likelihoods for model comparison (Bayes factors)
• Allows computation of posterior predictive distributions for model checking
• Handles high-dimensional parameter spaces common in hierarchical models

Algorithm structure

Proposal distribution

• Defines mechanism for generating candidate moves in parameter space
• Often symmetric (Gaussian) but can be asymmetric for improved efficiency
• Impacts mixing and convergence speed of the Markov chain
• Tuning proposal distribution crucial for algorithm performance

Acceptance probability

• Determines whether proposed move accepted or rejected
• Calculated using ratio of target densities and proposal densities
• Ensures detailed balance condition maintains correct stationary distribution
• Incorporates Metropolis-Hastings ratio $\alpha = \min\left(1, \frac{p(x')q(x|x')}{p(x)q(x'|x)}\right)$

Transition kernel

• Defines probability of moving from current state to proposed state
• Combines proposal distribution and acceptance probability
• Satisfies reversibility condition for Markov chain
• Ensures convergence to target distribution through repeated application

Implementation steps

Initial state selection

• Chooses starting point for Markov chain in parameter space
• Can impact convergence speed but not final equilibrium distribution
• Multiple chains with diverse starting points help assess convergence
• Often based on prior knowledge or maximum likelihood estimates

Proposal generation

• Draws candidate state from proposal distribution centered at current state
• Symmetric proposals (random walk) simplify acceptance probability calculation
• Adaptive proposals adjust scale based on acceptance rate during sampling
• Importance of balance between local and global exploration of parameter space

Acceptance-rejection mechanism

• Computes acceptance probability using Metropolis-Hastings ratio
• Generates uniform random number to compare against acceptance probability
• Accepts proposal if random number less than acceptance probability
• Retains current state if proposal rejected maintaining detailed balance

Convergence and diagnostics

Burn-in period

• Initial samples discarded to reduce impact of starting point
• Length determined by monitoring chain behavior and convergence diagnostics
• Ensures samples approximate target distribution more closely
• Can vary depending on complexity of target distribution and initial state choice

Autocorrelation assessment

• Measures dependence between successive samples in Markov chain
• High autocorrelation indicates slow mixing and potential convergence issues
• Effective sample size accounts for autocorrelation in uncertainty estimates
• Thinning (subsampling) reduces autocorrelation but may decrease efficiency

Gelman-Rubin statistic

• Compares within-chain and between-chain variances for multiple chains
• Values close to 1 indicate convergence to same target distribution
• Sensitive to differences in chain behavior and starting points
• Commonly used in conjunction with trace plots and other diagnostics

Advantages and limitations

Flexibility in complex models

• Handles intractable likelihood functions and complex prior distributions
• Applicable to high-dimensional parameter spaces and hierarchical models
• Enables sampling from non-standard distributions without closed-form solutions
• Facilitates exploration of multimodal and highly correlated posteriors

Computational efficiency considerations

• Can be computationally intensive for large datasets or complex models
• Efficiency depends on careful tuning of proposal distribution and algorithm parameters
• Parallel implementations and GPU acceleration improve performance
• Trade-off between exploration (acceptance rate) and exploitation (step size)

Variants and extensions

Metropolis algorithm

• Special case of Metropolis-Hastings with symmetric proposal distribution
• Simplifies acceptance probability calculation to ratio of target densities
• Often used with random walk proposals (Gaussian or uniform)
• Suitable for problems with well-behaved target distributions

Gibbs sampling

• Samples each parameter conditionally on others and current data
• Particularly efficient for models with conjugate prior-likelihood pairs
• Automatically accepts all proposals leading to faster convergence
• Can be combined with Metropolis-Hastings steps for non-conjugate components

Hamiltonian Monte Carlo

• Incorporates gradient information to propose more efficient moves
• Reduces random walk behavior and improves exploration of parameter space
• Particularly effective for high-dimensional and correlated posteriors
• Requires tuning of step size and number of steps in leapfrog integrator

Applications in Bayesian statistics

Parameter estimation

• Generates samples from posterior distribution of model parameters
• Enables calculation of point estimates (posterior mean, median) and credible intervals
• Handles complex joint posterior distributions with correlations and multimodality
• Facilitates inference on derived quantities and transformations of parameters

Model comparison

• Estimates marginal likelihoods for computing Bayes factors
• Implements reversible jump MCMC for comparing models with different dimensions
• Allows calculation of posterior model probabilities in Bayesian model averaging
• Enables assessment of predictive performance through cross-validation techniques

Posterior predictive checks

• Generates samples from posterior predictive distribution
• Assesses model fit by comparing replicated data to observed data
• Identifies systematic discrepancies between model and data
• Facilitates calculation of posterior predictive p-values for goodness-of-fit testing

Tuning and optimization

Proposal distribution selection

• Impacts efficiency and convergence rate of Metropolis-Hastings algorithm
• Gaussian proposals with adaptive covariance matrix common for continuous parameters
• Discrete proposals (uniform or customized) used for categorical variables
• Mixture proposals combine local and global moves for improved exploration

Acceptance rate targets

• Optimal acceptance rates depend on dimensionality and target distribution
• Generally aim for 20-50% acceptance rate for random walk Metropolis
• Higher acceptance rates (60-80%) preferred for Hamiltonian Monte Carlo
• Monitoring and adjusting proposal scale during sampling improves efficiency

Adaptive schemes

• Automatically tune proposal distribution during sampling
• Adjust scale and covariance of proposals based on acceptance history
• Implement diminishing adaptation to ensure ergodicity of Markov chain
• Balance exploration of parameter space with exploitation of high-density regions

Practical considerations

Dimensionality challenges

• Efficiency decreases as number of parameters increases (curse of dimensionality)
• Proposal distributions become less effective in high-dimensional spaces
• Requires careful parameterization and potentially dimension reduction techniques
• Specialized algorithms (Hamiltonian Monte Carlo) more robust to high dimensions

Multimodal distributions

• Standard Metropolis-Hastings can struggle with widely separated modes
• Parallel tempering introduces auxiliary chains at different temperatures
• Adaptive proposals or mixture proposals help jump between modes
• Importance of multiple chains and diverse starting points for thorough exploration

Parallel tempering

• Runs multiple chains at different "temperatures" to improve mixing
• Allows swaps between chains to escape local modes
• Gradually cools chains to focus on target distribution
• Particularly useful for complex, multimodal posteriors in phylogenetics and astrophysics

Software implementation

R packages

• MCMCpack provides functions for common Bayesian models
• rstan interfaces with Stan probabilistic programming language
• coda offers convergence diagnostics and output analysis tools
• bayesplot generates diagnostic plots and visualizations for MCMC output

Python libraries

• PyMC3 enables probabilistic programming and MCMC sampling
• emcee implements ensemble sampling algorithms for astrophysics
• TensorFlow Probability provides MCMC samplers with automatic differentiation
• arviz offers visualization and diagnostics for Bayesian inference

BUGS and JAGS

• Flexible environments for specifying and fitting Bayesian models
• Automatically choose appropriate MCMC algorithms based on model structure
• Widely used in epidemiology, ecology, and social sciences
• Integrate with R (R2jags, R2OpenBUGS) and other statistical software

Case studies and examples

Bayesian regression

• Estimates regression coefficients and variance parameters
• Incorporates prior information on coefficients and model structure
• Handles hierarchical models with varying intercepts and slopes
• Enables prediction and uncertainty quantification for new observations

Mixture models

• Infers number of components and component parameters simultaneously
• Implements label switching moves for exploring different clusterings
• Applies to density estimation, clustering, and classification problems
• Requires careful initialization and convergence assessment

Hierarchical models

• Pools information across groups while allowing for group-specific effects
• Estimates hyperparameters controlling distribution of group-level parameters
• Handles partial pooling between complete pooling and no pooling extremes
• Applied in meta-analysis, small area estimation, and longitudinal studies