In statistics, Markov chain Monte Carlo (MCMC) methods comprise a class of algorithms for sampling from a probability distribution. By constructing a Markov chain that has the desired distribution as its equilibrium distribution, one can obtain a sample of the desired distribution by recording states from the chain. The more steps that are included, the more closely the distribution of the sample matches the actual desired distribution. Various algorithms exist for constructing chains, including the Metropolis–Hastings algorithm.

Application domains

MCMC methods are primarily used for calculating numerical approximations of multi-dimensional integrals, for example in Bayesian statistics, computational physics,[1] computational biology[2] and computational linguistics.[3][4]

In Bayesian statistics, the recent development of MCMC methods has made it possible to compute large hierarchical models that require integrations over hundreds to thousands of unknown parameters.[5]

In rare event sampling, they are also used for generating samples that gradually populate the rare failure region.[citation needed]

General explanation

Convergence of the Metropolis–Hastings algorithm. Markov chain Monte Carlo attempts to approximate the blue distribution with the orange distribution.

Markov chain Monte Carlo methods create samples from a continuous random variable, with probability density proportional to a known function. These samples can be used to evaluate an integral over that variable, as its expected value or variance.

Practically, an ensemble of chains is generally developed, starting from a set of points arbitrarily chosen and sufficiently distant from each other. These chains are stochastic processes of "walkers" which move around randomly according to an algorithm that looks for places with a reasonably high contribution to the integral to move into next, assigning them higher probabilities.

Random walk Monte Carlo methods are a kind of random simulation or Monte Carlo method. However, whereas the random samples of the integrand used in a conventional Monte Carlo integration are statistically independent, those used in MCMC are autocorrelated. Correlations of samples introduces the need to use the Markov chain central limit theorem when estimating the error of mean values.

These algorithms create Markov chains such that they have an equilibrium distribution which is proportional to the function given.

Reducing correlation

While MCMC methods were created to address multi-dimensional problems better than generic Monte Carlo algorithms, when the number of dimensions rises they too tend to suffer the curse of dimensionality: regions of higher probability tend to stretch and get lost in an increasing volume of space that contributes little to the integral. One way to address this problem could be shortening the steps of the walker, so that it does not continuously try to exit the highest probability region, though this way the process would be highly autocorrelated and expensive (i.e. many steps would be required for an accurate result). More sophisticated methods such as Hamiltonian Monte Carlo and the Wang and Landau algorithm use various ways of reducing this autocorrelation, while managing to keep the process in the regions that give a higher contribution to the integral. These algorithms usually rely on a more complicated theory and are harder to implement, but they usually converge faster.


Random walk

Interacting particle methods

Interacting MCMC methodologies are a class of mean-field particle methods for obtaining random samples from a sequence of probability distributions with an increasing level of sampling complexity.[12] These probabilistic models include path space state models with increasing time horizon, posterior distributions w.r.t. sequence of partial observations, increasing constraint level sets for conditional distributions, decreasing temperature schedules associated with some Boltzmann–Gibbs distributions, and many others. In principle, any Markov chain Monte Carlo sampler can be turned into an interacting Markov chain Monte Carlo sampler. These interacting Markov chain Monte Carlo samplers can be interpreted as a way to run in parallel a sequence of Markov chain Monte Carlo samplers. For instance, interacting simulated annealing algorithms are based on independent Metropolis–Hastings moves interacting sequentially with a selection-resampling type mechanism. In contrast to traditional Markov chain Monte Carlo methods, the precision parameter of this class of interacting Markov chain Monte Carlo samplers is only related to the number of interacting Markov chain Monte Carlo samplers. These advanced particle methodologies belong to the class of Feynman–Kac particle models,[13][14] also called Sequential Monte Carlo or particle filter methods in Bayesian inference and signal processing communities.[15] Interacting Markov chain Monte Carlo methods can also be interpreted as a mutation-selection genetic particle algorithm with Markov chain Monte Carlo mutations.

Markov Chain quasi-Monte Carlo (MCQMC)

The advantage of low-discrepancy sequences in lieu of random numbers for simple independent Monte Carlo sampling is well known.[18] This procedure, known as Quasi-Monte Carlo method (QMC),[19] yields an integration error that decays at a superior rate to that obtained by IID sampling, by the Koksma–Hlawka inequality. Empirically it allows the reduction of both estimation error and convergence time by an order of magnitude.[citation needed] The Array–RQMC method combines randomized quasi–Monte Carlo and Markov chain simulation by simulating chains simultaneously in a way that the empirical distribution of the states at any given step is a better approximation of the true distribution of the chain than with ordinary MCMC.[20] In empirical experiments, the variance of the average of a function of the state sometimes converges at rate or even faster, instead of the Monte Carlo rate.[21]


Usually it is not hard to construct a Markov chain with the desired properties. The more difficult problem is to determine how many steps are needed to converge to the stationary distribution within an acceptable error.[22] A good chain will have rapid mixing: the stationary distribution is reached quickly starting from an arbitrary position. A standard empirical method to assess convergence is to run several independent simulated Markov chains and check that the ratio of inter-chain to intra-chain variances for all the parameters sampled is close to 1.[22][23]

Typically, Markov chain Monte Carlo sampling can only approximate the target distribution, as there is always some residual effect of the starting position. More sophisticated Markov chain Monte Carlo-based algorithms such as coupling from the past can produce exact samples, at the cost of additional computation and an unbounded (though finite in expectation) running time.

Many random walk Monte Carlo methods move around the equilibrium distribution in relatively small steps, with no tendency for the steps to proceed in the same direction. These methods are easy to implement and analyze, but unfortunately it can take a long time for the walker to explore all of the space. The walker will often double back and cover ground already covered.

Further consideration of convergence is at Markov chain central limit theorem. See [24] for a discussion of the theory related to convergence and stationarity of the Metropolis–Hastings algorithm.


Several software programs provide MCMC sampling capabilities, for example:

See also



  1. ^ Kasim, M.F.; Bott, A.F.A.; Tzeferacos, P.; Lamb, D.Q.; Gregori, G.; Vinko, S.M. (September 2019). "Retrieving fields from proton radiography without source profiles". Physical Review E. 100 (3): 033208. arXiv:1905.12934. Bibcode:2019PhRvE.100c3208K. doi:10.1103/PhysRevE.100.033208. PMID 31639953. S2CID 170078861.
  2. ^ Gupta, Ankur; Rawlings, James B. (April 2014). "Comparison of Parameter Estimation Methods in Stochastic Chemical Kinetic Models: Examples in Systems Biology". AIChE Journal. 60 (4): 1253–1268. doi:10.1002/aic.14409. PMC 4946376. PMID 27429455.
  3. ^ See Gill 2008.
  4. ^ See Robert & Casella 2004.
  5. ^ Banerjee, Sudipto; Carlin, Bradley P.; Gelfand, Alan P. (2014-09-12). Hierarchical Modeling and Analysis for Spatial Data (Second ed.). CRC Press. p. xix. ISBN 978-1-4398-1917-3.
  6. ^ Geman, Stuart; Geman, Donald (November 1984). "Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images". IEEE Transactions on Pattern Analysis and Machine Intelligence. PAMI-6 (6): 721–741. doi:10.1109/TPAMI.1984.4767596. ISSN 0162-8828. PMID 22499653. S2CID 5837272.
  7. ^ Gilks, W. R.; Wild, P. (1992-01-01). "Adaptive Rejection Sampling for Gibbs Sampling". Journal of the Royal Statistical Society. Series C (Applied Statistics). 41 (2): 337–348. doi:10.2307/2347565. JSTOR 2347565.
  8. ^ Gilks, W. R.; Best, N. G.; Tan, K. K. C. (1995-01-01). "Adaptive Rejection Metropolis Sampling within Gibbs Sampling". Journal of the Royal Statistical Society. Series C (Applied Statistics). 44 (4): 455–472. doi:10.2307/2986138. JSTOR 2986138.
  9. ^ Lee, Se Yoon (2021). "Gibbs sampler and coordinate ascent variational inference: A set-theoretical review". Communications in Statistics - Theory and Methods. 51 (6): 1–21. arXiv:2008.01006. doi:10.1080/03610926.2021.1921214. S2CID 220935477.
  10. ^ See Stramer 1999.
  11. ^ See Green 1995.
  12. ^ Del Moral, Pierre (2013). Mean field simulation for Monte Carlo integration. Chapman & Hall/CRC Press. p. 626.
  13. ^ Del Moral, Pierre (2004). Feynman–Kac formulae. Genealogical and interacting particle approximations. Springer. p. 575.
  14. ^ Del Moral, Pierre; Miclo, Laurent (2000). "Branching and Interacting Particle Systems Approximations of Feynman-Kac Formulae with Applications to Non-Linear Filtering". In Jacques Azéma; Michel Ledoux; Michel Émery; Marc Yor (eds.). Séminaire de Probabilités XXXIV (PDF). Lecture Notes in Mathematics. Vol. 1729. pp. 1–145. doi:10.1007/bfb0103798. ISBN 978-3-540-67314-9.
  15. ^ Del Moral, Pierre (2006). "Sequential Monte Carlo samplers". Journal of the Royal Statistical Society. Series B (Statistical Methodology). 68 (3): 411–436. arXiv:cond-mat/0212648. doi:10.1111/j.1467-9868.2006.00553.x. S2CID 12074789.
  16. ^ Chen, S.; Dick, Josef; Owen, Art B. (2011). "Consistency of Markov chain quasi-Monte Carlo on continuous state spaces". Annals of Statistics. 39 (2): 673–701. arXiv:1105.1896. doi:10.1214/10-AOS831.
  17. ^ Tribble, Seth D. (2007). Markov chain Monte Carlo algorithms using completely uniformly distributed driving sequences (Diss.). Stanford University. ProQuest 304808879.
  18. ^ Papageorgiou, Anargyros; Traub, J. F. (1996). "Beating Monte Carlo". Risk. 9 (6): 63–65.
  19. ^ Sobol, Ilya M (1998). "On quasi-monte carlo integrations". Mathematics and Computers in Simulation. 47 (2): 103–112. doi:10.1016/s0378-4754(98)00096-2.
  20. ^ L'Ecuyer, P.; Lécot, C.; Tuffin, B. (2008). "A Randomized Quasi-Monte Carlo Simulation Method for Markov Chains" (PDF). Operations Research. 56 (4): 958–975. doi:10.1287/opre.1080.0556.
  21. ^ L'Ecuyer, P.; Munger, D.; Lécot, C.; Tuffin, B. (2018). "Sorting Methods and Convergence Rates for Array-RQMC: Some Empirical Comparisons". Mathematics and Computers in Simulation. 143: 191–201. doi:10.1016/j.matcom.2016.07.010.
  22. ^ a b Gelman, A.; Rubin, D.B. (1992). "Inference from iterative simulation using multiple sequences (with discussion)" (PDF). Statistical Science. 7 (4): 457–511. Bibcode:1992StaSc...7..457G. doi:10.1214/ss/1177011136.
  23. ^ Cowles, M.K.; Carlin, B.P. (1996). "Markov chain Monte Carlo convergence diagnostics: a comparative review". Journal of the American Statistical Association. 91 (434): 883–904. CiteSeerX doi:10.1080/01621459.1996.10476956.
  24. ^ Hill, S. D.; Spall, J. C. (2019). "Stationarity and Convergence of the Metropolis-Hastings Algorithm: Insights into Theoretical Aspects". IEEE Control Systems Magazine. 39 (1): 56–67. doi:10.1109/MCS.2018.2876959. S2CID 58672766.
  25. ^ Foreman-Mackey, Daniel; Hogg, David W.; Lang, Dustin; Goodman, Jonathan (2013-11-25). "emcee: The MCMC Hammer". Publications of the Astronomical Society of the Pacific. 125 (925): 306–312. arXiv:1202.3665. Bibcode:2013PASP..125..306F. doi:10.1086/670067. S2CID 88518555.


Further reading