10. Equilibration, autocorrelation and uncertainty
Position in the course: Lesson 10 of 11. Complete the preceding derivation and use the explained exercises to check understanding.
1. A long file is not necessarily a long experiment
Trajectory frames are correlated and may all remain inside one metastable basin. Saving more frequently increases storage without necessarily increasing information. Before estimating an equilibrium quantity, define the target ensemble, an equilibration rule, the observable and a sampling interval. Stationarity is necessary for the formulas below; a drifting observable should not receive an equilibrium error bar merely by averaging it.
2. Derive the variance of an average
For \(M\) stationary observations \(A_n\) with variance \(\sigma_A^2\) and normalized autocorrelation \(\rho(k)\), expand every covariance in the average:
There are \(M-k\) pairs at lag \(k\) in each direction, which explains the finite-length weight. If correlations decay on scales much shorter than \(M\), define the statistical inefficiency and effective sample number:
Here the integrated autocorrelation time in saved-frame units is \(\tau_{int}=g/2\). Some texts use a convention that subtracts the half-lag contribution; report the convention. The autocorrelation tail is itself noisy and cannot be summed blindly to the full trajectory length.
3. A worked correlated sequence
For an analytic AR(1) model with \(\rho(k)=a^k\), \(|a|<1\), geometric summation gives \(g=(1+a)/(1-a)\). If \(a=0.8\), then \(g=9\). A record of \(9000\) values carries about \(1000\) independent-value equivalents, and the standard error is three times the independent-sample estimate. This calculation verifies the formula without generating random data. Negative correlations may give \(g<1\); the quantity depends on the observable, not just the integrator.
4. Blocking and independent replicas
Divide the stationary part into \(B\) blocks of length \(b\) and calculate each block mean. When \(b\) is large compared with correlation time, the block means become approximately independent and \(SE\) can be estimated from their scatter divided by \(\sqrt B\). Increase block length to seek a stable uncertainty plateau while retaining enough blocks to estimate its scatter. One giant block does not permit an error estimate.
For a derived quantity such as diffusion slope or a ratio, reanalyze whole blocks or independent replicas. Treating adjacent fit points as independent usually underestimates uncertainty. A block bootstrap must resample chunks long enough to retain relevant correlations. Independent starting states and replicas test trapping; blocks from one trapped trajectory cannot reveal an unvisited basin.
5. An error ledger
Keep statistical uncertainty separate from timestep bias, force tolerance, finite-cell effects and model uncertainty. A \(95\%\) confidence interval from one approximate ensemble does not include errors of a force field or functional. Choose convergence tolerances based on the target observable: density may equilibrate before a defect concentration or a conformational population. Define equilibration without repeatedly tuning the cutoff to obtain a preferred answer.
6. Exercises
How much longer must a stationary run be to halve its standard error if \(g\) is unchanged? Four times longer, because \(SE\propto M^{-1/2}\). Halving the output interval does not achieve this if neighboring frames merely become more correlated.
Two replicas have tiny within-replica error bars but incompatible means. What should be investigated? Different states, insufficient equilibration, metastability or inconsistent inputs. Combining them into a deceptively small error bar hides the evidence instead of resolving it.
7. References
- Chodera (2016): automated equilibration detection
- Rutgers computational physics: MD and Monte Carlo teaching material
Original analytic teaching illustration under the stated assumptions; no simulation data.
8. Related theory and practice
Quantum mechanics · Monte Carlo · Quantum Monte Carlo · CP2K · Quantum-Espresso