# MCMC autocorrelation convergence diagnostic

I use MCMC (Metropolis-Hastings) to sample posterior distributions of three parameters using a nonlinear least-squares objective function to calculate the likelihood of a parameter sets.

I want posteriors for about 500 different conditions, each time I calculate a posterior with MCMC.

Since it is hard to check all 500 runs by hand for convergence, I would like to have a simple rule based on autocorrelation which indicates for each run if the chain converged or not, checking all conditions individually would be very cumbersome.

Each chain has length 40000 (burn-in already removed) and I calculated autocorrelations for lags 1 to 10000.

I then calculated the autocorrelation time as the time constant for the autocorrelation series (the first lag when the autocorrelation drops below 1/exp(1)).

I get autocorrelation times between 16 and about 7000.

I suppose that with an autocorrelation time of 16 I can be sure that the Markov chain converged, but is there maybe some "rule of thumb" which tells me, given a chain of length 40000, which autocorrelation time is still "acceptable"?

I searched and looked for instance here, however, the approaches I found require a visual investigation.

• I believe what most people do (if the algorithms are fast enough and time is no issue) is to take the highest boundry, that is, judge the worst-behaving draw (in terms of autocorrelation) by visual inspection, multiply the resulting burn-in by 2, and apply that to all chains. Commented Jul 5, 2013 at 14:06
• Does this mean I calculate the new burn-in time using the autocorrelation and then calculate my chains again? Time is an issue for me, so I would like to assess convergence from the chains that I have now... Commented Jul 5, 2013 at 15:00

The Geweke's convergence diagnostic works with a single chain and provides a z-score. There is also the Heidelberger and Welch's convergence diagnostic, which also works with a single chain and gives a p-value of the null hypothesis that the sampled values come from a stationary distribution. If you run multiple chains, there's Gelman and Rubin's convergence diagnostic. If you use R, all of these are implemented in the coda package.
• Thanks for the quick answer, I looked into Geweke's convergence diagnostic in coda and calculated the z-scores. Then I calculated p-values using 2*pnorm(-abs(z)) Commented Jul 5, 2013 at 14:53