8
$\begingroup$

There is a lot of literature out there about Markov chain Monte Carlo (MCMC) convergence diagnostics, including the most popular Gelman-Rubin diagnostic. However, all of these assess the convergence of the Markov chain, and thus address the question of burn-in.

Once I have figured out burn-in, how should I decide how many MCMC samples are enough to continue with my estimation process? Most of the papers using MCMC mention that they ran the Markov chain for some $n$ iterations, but do not say anything about why/how they chose that number, $n$.

In addition, one desired sample size cannot be the answer for all samplers, since the correlation in the Markov chain varies greatly from problem to problem. So is there a rule out there to find out the number of samples required?

|cite|improve this question
$\endgroup$
4
$\begingroup$

How many samples (post burn-in) that you need depends on what you are trying to do with those samples and how your chain mixes.

Typically we are interested in posterior expectations (or quantiles) and we approximate these expectations by averages of our posterior samples, i.e. $$ E[h(\theta)|y] \approx \frac{1}{M} \sum_{m=1}^M h(\theta^{(m)}) = E_M $$ where $\theta^{(m)}$ are samples from your MCMC. By the Law of Large Numbers the MCMC estimate converges almost surely to the desired expectation.

But to address the question of how many samples we need to be assured that we are close enough to the desired expectation, we need a Central Limit Theorem (CLT) result, i.e. something like $$ \frac{E_M -E[h(\theta)|y]}{\sqrt{M}} \stackrel{d}{\to} N(0,v_h^2) $$ With this CLT we could then make probabilitic statements such as "there is 95% probability that $E[h(\theta)|y]$ is between $E_M \pm 1.96 v_h$." The two issues here are

  1. Does the CLT apply?
  2. How can we estimate $v_h^2$.

We have some results on when the CLT applies, e.g. discete state chains, reversible chains, geometrically ergodic chains. See Robert and Casella (2nd ed) section 6.7.2 for some results in this direction. Unfortunately the vast majority of Markov chains that are out there have no proof that CLT exists.

If a CLT does exist, we still need to estimate the variance in the CLT. One way to estimate this variance involves breaking the chain up into blocks, see Gong and Flegal and references therein. The method has been implemented in the R package mcmcse with the functions mcse and mcse.q to estimate the variance for expectations and quantiles.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ This sounds reasonable, and I am familiar with the literature here. How often is this argument actually used in practice? $\endgroup$ – Greenparker Mar 23 '16 at 18:32
2
$\begingroup$

John Kruschke in Doing Bayesian Data Analysis recommends that for parameters of interest, MCMC chains should be run until their effective sample size is at least 10,000. Although no simulations are presented, I believe his rationale is that ESS > 10,000 ensures numerically stable estimates. However, I have seen that Andrew Gelman and other Stan developers recommend less (e.g. 2000 - 3000 is fine; no exact link, unfortunately, but see discussions on the Stan Google users group). Moreover, for complex models, running chains long enough for an ESS > 10,000 can be arduous!

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks. Can you send me where he says that in his material? It'll take a long time to skim through the webpage. Also, my answer [here] talks about determining lower bound for ESS. $\endgroup$ – Greenparker Jul 21 '16 at 14:29
  • $\begingroup$ Sorry, I realized I didn't put the link. Here it is. $\endgroup$ – Greenparker Jul 21 '16 at 14:37
  • 1
    $\begingroup$ Sorry, I should have been more specific. Kruschke mentions it briefly in his blog post here doingbayesiandataanalysis.blogspot.co.uk and it is in Chapter 7 of his book, 'Markov Chain Monte Carlo', page 184 of the 2nd edition: books.google.co.uk/…. $\endgroup$ – user3237820 Jul 21 '16 at 15:26
1
$\begingroup$

This is indeed one of the big downsides of MCMC simulations, since there is no general and apriori estimate on the number of samples. I think a good literature here is "Some things we've learned (about MCMC)" by Persi Diaconis which deals with a lot of subtleties of MCMC simulations which might indicates that there is no easy answer to this question.

In general, making good estimates on how long to run the chain requires a good understanding of the mixing time of the Markov chain, which heavily depends on properties of the underlying graph. There are many "burn-in-free"-methods out there to bound the mixing time from above, but all of these methods have in common that they need a deeper understanding of the underlying Markov chain and the involved constants are typically hard to compute. See for example "Conductance and Rapidly Mixing Markov Chains" by King, "Path Coupling: a Technique for Proving Rapid Mixing in Markov Chains" by Bubley et al., or "Nash inequalities for finite Markov chains" by Diaconis et al.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Agreed. But in practice the mixing time of samplers is not always studied in such detail to be address this question. In addition, studying the mixing time requires considerable background in Markov chain theory, something that most end users of MCMC might not be familiar with. There are no even heuristics out there (like the diagnostics)? $\endgroup$ – Greenparker Mar 23 '16 at 17:54
  • $\begingroup$ The only thing I can think of is to numerically estimate the second largest eigenvalue of the transition matrix and to derive a bound on the mixing time from that. You may have a look at the PhD thesis of Kranthi Kumar Gade. $\endgroup$ – Tobias Windisch Mar 23 '16 at 18:01
  • $\begingroup$ What if I am working with a general state space Markov chain, not a finite state space? I am guessing that is not possible then, but will take a look at it. $\endgroup$ – Greenparker Mar 23 '16 at 18:10
  • $\begingroup$ You are right. His method only works for finite state spaces and discrete-time Markov chains. $\endgroup$ – Tobias Windisch Mar 23 '16 at 18:13

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.