# Convergence tests in Markov Chain Monte Carlo

For a relatively simple Markov chain Monte Carlo process, such as using Metropolis to find calculate thermal averages for an Ising model, how is it possible to determine whether quantities have converged?

If one knows the autocorrelation time, this seems relatively simple. You just run it for a sufficiently high multiple of this time and Bob's your uncle.

If you don't know the autocorrelation time, it would seem more complex. No matter how hard you try to ensure that your averages have converged, there's always the possibility that it is stuck in a local minima. This is especially troublesome if the autocorrelation time scales polynomially or even exponentially with the system size.

So how can you ensure that the quantities really have converged? What convergence tests are not fooled by local minima? I read here recently about the method of logarithmic binning. Could that do the job?

-

## migrated from physics.stackexchange.comJun 8 '13 at 22:09

This question came from our site for active researchers, academics and students of physics.

I think there are almost by definition no such tests. Suppose you use a 'bad' (non-cluster) algorithm on the Ising model in the broken phase, then indeed you can simulate for days getting wrong estimates for your observables. However, your data will be perfectly consistent and no data analysis/binning/whatever will tell you that you're wrong, unless you extend the runtime of your simulation. – Vibert Jun 8 '13 at 20:49
Yes, it's a difficult problem. Nevertheless, in simple enough situations, such as the one mentioned here, it is possible to implement a perfect simulation algorithm, i.e., one that returns (after a random, but almost surely finite, time) a configuration sampled exactly from the Gibbs measure. Relevant keywords are: Propp-Wilson algorithm and coupling from the past. – Yvan Velenik Jun 8 '13 at 21:53

It's very hard to prove convergence for an ensemble in statistical mechanics, for exactly the reasons you've suggested—you may be stuck moving around a local minimum, instead of sampling the complete phase space. About the best that you can do is to demonstrate likely convergence of your simulation.

However, you can do calculations to demonstrate autocorrelation of various properties. For instance, you can measure the autocorrelation of the energy function: $\left<E(t)E(t + \delta t) \right>$. Alternatively, you can use the multiscale error correlation method of Flyvbjerg and Petersen to determine when your data samples are effectively decorrelated from one another, and use that as a means of determining whether you have enough data (and to estimate the statistical error in the mean of your data samples).

-
How can such a correlation method 'likely' prove convergence? How does it know about other local minima your data are not seeing? – Vibert Jun 10 '13 at 13:15
For starters, if you're stuck in a deep minimum, you're unlikely to get away from it. Hence your decorrelation time, as measured by the algorithm, would be quite long, relative to the length of your simulation. That would mean you've sampled relatively few independent states, which would violate the general principles behind molecular simulations (which states that large numbers of data points, generated either by time trajectories or random sampling of phase space will generate a good estimate of the "true" observable). – aeismail Jun 10 '13 at 16:17

I can't comment yet, so I would be grateful if someone could turn this into a comment to aeismail's answer. Although it's quite long.

I don't have access to the above mentioned Flyvbjerg and Petersen paper and thus can't really judge on what follows, but in the really nice paper "Quantifying uncertainty and sampling quality in biomolecular simulations" by Grossfield and Zuckerman, the authors comment: "The approach can be described simply (although it is not easily understood from the original reference)." In my opinion, Grossfield and Zuckerman give a really nice and simple explanation of this procedure which they call "block averaging" and how to obtain estimates of correlation times with it.

Another paper which might be interesting is Automated Sampling Assessment for Molecular Simulations Using the Effective Sample Size by Zhang, Bhatt and Zuckerman. In this article, they assess sampling quality and convergence via an Effective Sample Size (ESS) estimated from the variances in the populations of physical states. They argue that state populations are the most critical slow observables and thus are a good observable to estimate correlation times from. But they also put the fundamental problem in a kind of disclaimer, stating that they answer the question: "What is the statistical quality of the sampling based on the configurational states visited in a given set of simulations?" and they of course are aware that "without prior knowledge or assumptions about a landscape, it would appear impossible to know whether every important state has been visited in a given simulation."

-