How to interpret autocorrelation plot in MCMC

Question

I am getting familiar with Bayesian statistics by reading the book Doing Bayesian Data Analysis, by John K. Kruschke also known as the "puppy book". In chapter 9, hierarchical models are introduced with this simple example: \begin{align} y_{ji} &\sim {\rm Bernoulli}(\theta_j) \\ \theta_j &\sim {\rm Beta}(\mu\kappa, (1-\mu)\kappa) \\ \mu &\sim {\rm Beta}(A_\mu, B_\mu) \\ \kappa &\sim {\rm Gamma}(S_\kappa, R_\kappa) \end{align} and the Bernoulli observations are 3 coins, each 10 flips. One shows 9 heads, the other 5 heads and the other 1 head.

I have used pymc to infer the hyperparamteres.

with pm.Model() as model:
# define the     
    mu = pm.Beta('mu', 2, 2)
    kappa = pm.Gamma('kappa', 1, 0.1)
    # define the prior
    theta = pm.Beta('theta', mu * kappa, (1 - mu) * kappa, shape=len(N))
    # define the likelihood
    y = pm.Bernoulli('y', p=theta[coin], observed=y)

    # Generate a MCMC chain
    step = pm.Metropolis()
    trace = pm.sample(5000, step, progressbar=True)
    trace = pm.sample(5000, step, progressbar=True)


burnin = 2000  # posterior samples to discard
thin = 10  # thinning 
pm.autocorrplot(trace[burnin::thin], vars =[mu, kappa])

My question is regarding the autocorrelation. How shall I interpret autocorrelation? Could you please help me to interpret the autocorrelation plot?

It says as the samples gets further from each other the correlation between them reduces. right? Can we use this to plot to find the optimum thinning? Does the thinning affect the posterior samples? after all, what is the use of this plot?

Answer 1

First of all: if memory and computational time for handling the MCMC output are not limiting, thinning is never "optimal". At an equal number of MCMC iterations, thinning of the chain always leads (on average) to a loss precision of the MCMC approximation.

Routinely thinning based on autocorrelation or any other diagnostics is therefore not advisable. See Link, W. A. & Eaton, M. J. (2012) On thinning of chains in MCMC. Methods in Ecology and Evolution, 3, 112-115.

In every-day practice, however, there is the common case that you have to work with a model for which the sampler doesn't mix very well (high autocorrelation). In this case

1) Close chain elements are very similar, meaning that throwing one away doesn't loose a lot of information (that is what the autocorrelation plot shows)

2) You need a lot of repetitions to get convergence, meaning that you get very large chains if you don't thin. Because of that, working with the full chain can be very slow, costs a lot of storage, or even lead to memory problems when monitoring a lot of variables.

3) Additionally, I have the feeling (but that one I have never tested systematically) that thinning makes JAGS a bit faster as well, so might be able to get a few more iterations in the same time.

So, my point is: the autocorrelation plot gives you a rough estimate about how much information you are loosing through thinning (note though that this is an average over the whole posterior, loss may be higher in particular regions).

Whether this price is worth paying depends on what you gain by thinning in terms saving computing resources and time later. If MCMC iterations are cheap, you can always compensate the loss of thinning by running a few more iterations as well.