MCMC; Can we be sure that we have a ''pure'' and ''large enough'' sample from the posterior? How can it work if we are not?

Question

Referring to this thread: How would you explain Markov Chain Monte Carlo (MCMC) to a layperson?.

I can see that it is a combination of Markov Chains and Monte Carlo: a Markov chain is created with the posterior as invariant limiting distribution and then Monte Carlo draws (dependent) are made from the limiting distribution (=our posterior).

Let's say (I know that I am simplifying here) that after $L$ steps we are at the limiting distribution $\Pi$ (*).

The Markov chain being a sequence of random variables, I get a sequence $X_1, X_2, \dots , X_L, \Pi, \Pi, \Pi, \dots \Pi$, where $X_i$ is a random variable and $\Pi$ is the limiting ''random variable'' from which we wish to sample.

The MCMC starts from an initial value, i.e. $X_1$ is a random variable with all mass at that one value $x_1$. If I use capital letters for random variables and small letters for realizations of a random variable, then the MCMC gives me a sequence $x_1,x_2,x_3, \dots x_L, \pi_1, \pi_2, \pi_3, ....\pi_n$. So the length of the MCMC chain is L+n.

[[*Note: the capital letters are random variables (i.e. a whole bunch of outcomes) and the small $x$ are outcomes, i.e. one particular value. * ]]

Obviously, only the $\pi_i$ belong to my ''posterior'' and for approximating the posterior ''well'' the value of $n$ should be ''large enough''.

If I summarise this then I have an MCMC chain $x_1,x_2,x_3, \dots x_L, \pi_1, \pi_2, \pi_3, ....\pi_n$ of length $N=L+n$, only $\pi_1,\pi_2,\dots, \pi_n$ are relevant for my posterior approximation, and $n$ should be large enough.

If I do include some of the $x_i$ (i.e. realizations before the invariant distribution is reached) in the computation of the approximation of the posterior, then it will be ''noisy''.

I know the length of the MCMC chain $N=L+n$, but without knowledge of the $L$, i.e. the step where I am sure to sample from the limiting distribution, I can not be sure that I did not include noise, nor can I be sure about $n=N-L$, the size of my sample from the limiting distribution, in particular, I can not be sure whether it is ''large enough''.

So, as far as I understood, this value of $L$ is of critical importance for the quality of approximation of the posterior (exclusion of noise and a large sample from it).

Are there any ways to find a reasonable estimate for $L$ when I apply MCMC ?

(*) I think that, in general, $L$ will depend on the initial value $x_1$.

Answer 1

TL DR; You can't estimate $L$ since $L = \infty$. Thus, the simplifying assumption can never truly be possible. (There maybe some cases where it is, but not in the general world of MCMC). You can however decide what $N$ will make the early bias small.

---

Essentially, your question boils down to "how can we estimate burn-in time?". Burn-in is the act of throwing away the beginning samples because the Markov chain has not converged. There are many MCMC diagnostics that help you estimate "burn-in" time, you can see a review of them here.

There are two schools of through regarding burn-in; the popular one is to use one of those diagnostics to decide what $L$ is, and throw away the $L$ samples, and the second school of through it, the first $L$ samples shouldn't matter, so don't worry about them. Charlie Geyer has a rant about this that I agree with.

Now, I turn to the more technical details of your question.

A simplifying assumption you make in your question is that eventually, (after $L$ steps), the sampler will start drawing from the limiting distribution. So your samples after $L$ steps are pure draws, albeit correlated. This is untrue. Strictly speaking, $L$ is $\infty$. The Markov chain never truly converges to the limiting distribution in finite time. So estimating $L$ is almost pointless.

A different way of posing this question is: what is $L$ such that after $L$ steps, the Markov chain is "close enough" to the limiting distribution. This is the question most diagnostics try to answer. It is increasingly agreed upon that the diagnostics above are generally extremely liberal, and can diagnose "convergence" much before it should have. Here is a paper that demonstrates some of the weaknesses of diagnostics.

What the above asks the users to do instead is don't worry about $L$, worry about $N$. Generally, users are not interested in the full posterior distribution, but in a specific quantity. Often this quantity is the mean of the posterior, or any other function that can be written down as an expectation. This is where the "Monte Carlo" part of MCMC comes in, since Monte Carlo indicates estimating an integral with summation. So if $X_1, X_2, X_3, \dots, X_N$ is your Markov chain (notice how I am ignoring $L$, since $L$ is $\infty$), and we want to estimate the posterior mean ($\theta$), then $$ \bar{\theta}_N = \dfrac{1}{N} \sum_{i=1}^{N}X_i. $$

The idea is that if $N$ is large enough, then the initial bias of the sample will be insignificant. Of course if the starting value was pathetically far away from high probability space of the limiting distribution, a user can eye-ball and throw away the first couple of samples. This is different from estimating $L$, since it is not an estimation, but an educated disregard for clearly corrupted samples.

Now the question of course is: how large should $N$ be? The answer should depend on how well do we want to estimate $\theta$. If we want a great estimate, then we want more samples, if an ok estimate suffices, then we might be fine with a smaller sample. This is also exactly what happens in standard statistical problems.

The way we quantify "goodness" of an estimate, is to think, "what can we say about $(\bar{\theta}_N - \theta)$, the Monte Carlo error? Under reasonable conditions, there is in fact a Markov chain CLT that says as $N \to \infty$, for any initial distribution $$\sqrt{N}(\bar{\theta}_N - \theta) \overset{d}{\to} N_p(0, \Sigma), $$

where $\theta \in \mathbb{R}^p$ and $\Sigma$ is the asymptotic covariance matrix. The key here is that the result is true for any initial distribution.

When $\Sigma/N$ is small, we know that the estimator is good. This paper presents this idea of stopping, and my answer here summarizes their method. The results in their paper are also regardless of the initial distribution of the process.