1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – mzuba Jul 5 '13 at 14:06
  • $\begingroup$ 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... $\endgroup$ – user1981275 Jul 5 '13 at 15:00
3
$\begingroup$

Have you looked at the standard tests for convergence?

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.

$\endgroup$
  • $\begingroup$ 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)) $\endgroup$ – user1981275 Jul 5 '13 at 14:53
  • $\begingroup$ I calculated the pvalues using Heidelberger and Welch as well and they differ quite substantially from the pvalues I calculated from the z-score from Geweke's diagnostic. Also, I see that if the autocorrelation time is very low (<100), the p-values are very high. I must still be missing something... $\endgroup$ – user1981275 Jul 5 '13 at 15:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.