Autocorrelation and Effective Sample Size

Autocorrelation changes the uncertainty of an average even when every stored configuration has the correct stationary distribution. Effective sample size translates that variance inflation into the number of independent observations with the same precision. It belongs to an estimator, not to a Markov chain in isolation.

Required background. MCMC stationarity and mixing distinguishes equilibration from production sampling.

Variance of a correlated average

Let $A_t$ be a real stationary observable with finite nonzero variance, $C(t)=\operatorname{Cov}(A_0,A_t)$, and $\rho(t)=C(t)/C(0)$. Expanding the variance of a sample mean gives the exact finite-$N$ identity

$$\operatorname{Var}(\bar A) =\frac{C(0)}N \left[1+2\sum_{t=1}^{N-1} \left(1-\frac tN\right)\rho(t)\right].$$

When correlations are summable and $N$ is much longer than their decay, define

$$\tau_{\rm int}=\frac12+\sum_{t=1}^{\infty}\rho(t), \qquad N_{\rm eff}=\frac{N}{2\tau_{\rm int}}.$$

For $\tau_{\rm int}>0$, $\operatorname{Var}(\bar A)\approx C(0)/N_{\rm eff}$. If the summed covariance is zero, the leading $1/N$ variance vanishes and the finite-$N$ expression must be used instead of interpreting an infinite effective count literally. This convention gives $\tau_{\rm int}=1/2$ for independent observations. Some literature defines twice this quantity as the integrated time; always state the convention.

An autoregressive benchmark

For $A_{t+1}=rA_t+\eta_t$ with $|r|<1$ and independent zero-mean innovations chosen for stationarity, $\rho(t)=r^t$. Therefore

$$2\tau_{\rm int}=\frac{1+r}{1-r}.$$

At $r=0.9$, the variance factor is $19$. Ten thousand stored observations provide approximately $526$ independent-observation equivalents, and the naive standard error is too small by $\sqrt{19}\approx4.36$. This is a useful synthetic test for both autocorrelation and blocking routines.

Negative $r$ yields anticorrelation and a factor below one. Algorithms that force every estimated correlation positive will give a conservative but distorted answer for this example.

Truncating the noisy correlation sum

An empirical correlation at long lag has large noise because few pairs remain. Summing it out to the full trajectory length does not improve the estimator. Use a justified window that extends past the visible correlation decay while remaining short relative to the run; compare nearby windows and independent replicates.

Long tails or multiple relaxation scales complicate this choice. An energy observable may decay rapidly while a phase-switching variable scarcely moves. A short estimated $\tau_{\rm int}$ in one metastable region does not diagnose equilibration across regions.

For nonlinear or reweighted observables, analyze the relevant combined fluctuation, such as $w(A-\mu)$ for a weighted ratio, rather than substituting the autocorrelation time of energy by habit.

Blocking as a complementary check

Divide a production sequence into $B$ contiguous blocks of length $L$, and compute each block mean $a_b$. If $L$ exceeds the relevant correlation time and enough blocks remain,

$$\widehat{\operatorname{SE}}(\bar A) =\sqrt{\frac{1}{B(B-1)} \sum_{b=1}^B(a_b-\bar A)^2}.$$

Increase $L$ and look for a stable error plateau. At very large $L$, too few blocks remain to estimate uncertainty reliably; a noisy plateau based on two or three blocks is not persuasive. Check independent chains as well.

Thinning reduces storage but usually throws away information when observations are cheap. It can be useful when measuring an observable is expensive, but the appropriate comparison is total computational cost at fixed precision.

Reporting limits

State update units, production length, discarded equilibration segment, block sizes or correlation windows, number of independent chains, and the observable-specific estimated uncertainty. If the run is not much longer than the slowest measured time, report that limitation rather than printing extra digits.

Exercises

A run records one value after every ten sweeps. Its estimated $\tau_{\rm int}$ is four in recorded-sample units. What is it in sweep units?

Solution

Four recorded-sample intervals correspond to forty sweeps. This unit conversion does not recover the integrated autocorrelation time of a series measured every sweep: the unrecorded lags are unknown. For the stored mean, use $N_{\rm eff}=N/(2\times4)$ with $N$ the stored count. Do not multiply the time by ten while leaving $N$ in recorded-sample units; that would mix units and overinflate the variance.

References

• Flyvbjerg, H., and H. G. Petersen, “Error estimates on averages of correlated data,” The Journal of Chemical Physics 91, 461–466 (1989), doi:10.1063/1.457480 — blocking correlated measurements.
• Hastings, W. K., “Monte Carlo sampling methods using Markov chains and their applications,” Biometrika 57, 97–109 (1970), doi:10.1093/biomet/57.1.97 — stationary sampling and proposal correction.