# Good method for correlated samples and estimating autocorrelation times

I'm working on a Monte Carlo project similar to the Ising model. I've found many examples on which I've based my code.

From some papers I read on binning analysis, the errors after each binning step are supposed to converge. Mine ended up oscillating after some binning step. And so I'm getting negative auto correlation times.

I was hoping someone could either verify my procedure is correct, or explain a good procedure for dealing with correlated sampling.

The algorithm I am using is as follows:

• I sample after each time step.

• Then recursively (the function in the code is binning()),

• I create a list of lists.

• Each list in the list is a collection of averages of successive samples by binning successive pairs (not crossing between pairs).

• I calculate the standard deviation based on a paper at each step and that gives a sequence of errors at each binning step.

Part 4D of this paper were my references. I assumed the last error would be a good estimate for the limit, giving me negative correlation as it doesn't converge, my problem.

I'm not sure if I am doing this correctly, and if I am, what could be going wrong giving somewhat oscillating error for each bin?

• Can you explain you procedure a bit more? – nicoguaro Feb 16 at 16:57
• @nicoguaro Sure. I sample after each time step. Then recursively (the function in the code is binning()), I create a list of lists. Each list in the list is a collection of averages of successive samples by binning successive pairs (not crossing between pairs). I calculate the standard deviation based on a paper at each step and that gives a sequence of errors at each binning step. Part 4D of this paper were my references: arxiv.org/pdf/0906.0943.pdf. I assumed the last error would be a good estimate for the limit, giving me negative correlation as it doesn't converge, my problem – Daniel Sela Feb 16 at 18:07
• You could add those details in the question ;). – nicoguaro Feb 16 at 19:38
• I mention it because the purpose of this site is not reviewing code, so, if you add an explanation it is easier that the question is on-topic. – nicoguaro Feb 16 at 19:39
• @nicoguaro Sorry, this is my first time posting here. I added the explanation. If there's anything else you think I should add, please let me know. – Daniel Sela Feb 16 at 19:57

I recently went through a similar experience, so you can take this suggestion as coming from a fellow learner rather than an expert. I was implementing an MCMC-type algorithm, and I wanted to know what the effective sample size was. So I tried to calculate the empirical autocorrelation $$\hat\rho_k$$ of the samples, but the sum

$$\tau = 1 + 2\sum_{k = 1}^N\hat\rho_k$$

was close to zero or negative. This would imply a really large effective sample size (ESS) if you took that estimate seriously, which was obviously absurd. The reason this works poorly is that it's difficult to get a very good estimate for the autocorrelation function at large lags from the samples themselves.

I found an arXiv preprint that compared several different methods for estimating the ESS and/or correlation time for functions of the samples. The easiest method was a binning approach similar to but I think simpler than the one in the paper you reference. First, let $$f$$ be some function of the samples that you're trying to estimate, let

$$\hat\mu = \frac{1}{N}\sum_{k = 1}^Nf(x_k)$$

be the sample mean of $$f$$, and

$$\hat s^2 = \frac{1}{N}\sum_{k = 1}^N(f(x_k) - \hat\mu)^2$$

the sample variance. Now divide the $$N$$ samples into batches of size $$M$$. (Usually $$M \propto N^{1/3}$$, but all you really need is that both $$M$$ and $$N / M$$ to go to infinity in the limit of larger sample sizes.) Now let

$$\hat\mu_m = \frac{1}{M}\sum_{k = M\cdot m + 1}^{M\cdot(m + 1)}f(x_k)$$

be the sample mean of batch $$m$$, and let

$$\hat s_M^2 = \frac{M}{N}\sum_{m = 1}^{N / M}(\hat\mu_m - \hat\mu)^2$$

be the variance between the batches. Then the correlation time is approximately

$$\tau = M\hat s_M^2 / \hat s^2.$$

I found this approach to be pretty foolproof to implement and the results are guaranteed to be positive, whereas the empirical correlation function can be poorly sampled and thus might sum to all sorts of crazy things.

Section IV.D of the paper you linked is similar to the method described in section 2.3 of the Thompson preprint, the initial positive sequence or IPS estimator. The idea is that the sequence

$$\textrm{IPS}_k = \rho_{2k} + \rho_{2k + 1}$$

is always positive for the true, ideal autocorrelation function $$\rho_k$$ of a reversible Markov chain. The idea of the IPS estimator is to take the correlation time to be

$$\tau \approx 1 + 2\sum_{k = 1}^{N^*}\textrm{IPS}_k$$

but this sum is truncated at the first index $$N^*$$ where $$\textrm{IPS}_{N^* + 1} < 0$$. I also tried this approach for my MCMC simulation and found that it gave about the same estimate of the correlation time and thus the effective sample size.

Finally, you might have better luck asking this kind of question on the Statistics stack exchange as it's a little more their bailiwick.

• Thank you for the long explanation and the paper you linked. Is there any concern about taking the batch size to be too small or too large. Or should this be done for multiple batches? Also, how does the autocorrelation time estimate the effective sample size? I will also try on the Statistics page, thanks for all the help, I've been stick on this for a while before I decided to ask here. – Daniel Sela Feb 16 at 22:21
• The effective sample size is equal to the number $N$ of samples divided by the correlation time $\tau$. When the samples are all uncorrelated, $\tau = 1$ and so the effective sample size is equal to $N$. Usually $\tau$ is much greater than 1. (There's a trick called antithetic Markov updates that makes successive estimates negatively correlated, giving an effective sample sizes that's greater than $N$. But that's deep hackery, above my pay grade.) – Daniel Shapero Feb 16 at 23:29
• That makes sense. I'd assume you want to choose N large enough so that the effective sample size is sufficiently large and 𝜏 sufficiently small. For the batching process you mentioned, is there any concern about taking M to be too small or too large? And lastly, is the actual error calculated as in equation 18 of the reference I attached using the 𝜏 as you have it from the batching? Sorry for all the questions. This is my first MCMC, so I am a bit inexperienced in this part of the estimation in the simulation. – Daniel Sela Feb 16 at 23:57
• No worries -- I found Practical MCMC by Geyer to be a great reference if you want to read more. I think as long as both the batch size $M$ and the number of batches $N / M$ go to infinity in the limit of large $N$, you're fine. In principle you could do $M = N^\alpha$ for any $\alpha$ between 0 and 1. You want a big enough batch size that you can get a good enough estimate of the batch mean, and enough batches to get a good estimate of the inter-batch variance. Both should probably be over 30 at least. – Daniel Shapero Feb 17 at 3:41
• I've implemented this and have gotten some times less than 1 and others above 1 (at different temperatures). Is this problematic because than the effective sampling size could be greater than the sampling size? – Daniel Sela Feb 17 at 4:51