Lab: error bars by hand

The setup

An average from a simulation is itself a random number: run again with other random velocities and you get a slightly different value. The error bar says how much. Two facts make MD special:

Notation: s is the standard deviation of the samples, "error" means the standard error of the mean (one σ, about 68 % confidence).

Step 1

The textbook error bar

Why: Start with the case where the formula is right: independent measurements. The spread of the samples tells how much one measurement varies; the error of the mean is smaller by √N because random deviations partly cancel.

What the numbers say: One reading varies by ±0.89, the mean of ten by ±0.28. To halve the error you need four times more independent samples: errors fall only as 1/√N. This is the formula most analysis scripts use by default, and for MD data it gives error bars that are far too small, as the next steps show.

Step 2

Correlated data: the statistical inefficiency

Why: In MD, a frame looks much like the one before it. Such samples carry less information than independent ones. The autocorrelation ρ(k) measures how similar values k frames apart are; adding it up tells how many frames count as one independent sample.

What the numbers say: Sixteen correlated numbers are worth about six independent ones, and the honest error bar is 1.66 times the naive one. For real MD, where frames may be saved every 10 fs while the system needs picoseconds to forget, g can be in the hundreds. The limitation: ρ(k) from a short series is noisy (here it swings to −0.56 at lag 3), so the sum must be cut off sensibly (at the first negative value is a common rule), and the result is only an estimate.

Step 3

Block averaging by hand

Why: Block averaging avoids computing ρ(k) at all. Cut the series into blocks, average each block, and treat the block means as samples. Once the blocks are longer than the memory of the data, the block means are independent and the naive formula becomes correct for them.

What the numbers say: The error grows with block size (0.095 → 0.130 → 0.147) and stops growing once blocks are longer than the correlation: the plateau, 0.147, agrees with the g-corrected 0.158. Beyond that, too few blocks remain (two blocks of 8 give 0.125, itself very uncertain). In practice: increase the block size until the error levels off, and read the plateau while there are still at least 20–50 blocks. Squaring the ratio plateau/naive gives g directly: (0.147/0.095)² = 2.4.

Step 4

A test with a known answer

Why: To trust the method, apply it where the exact answer is known: x(t+1) = φ x(t) + noise. Each value keeps a fraction φ of the previous one, the correlation is exactly ρ(k) = φᵏ, and g = (1 + φ)/(1 − φ).

What the numbers say: The block curve rises from the naive 0.0072 to a plateau of about 0.031 at blocks of 64–512, exactly the exact error 0.0308 and the real scatter of 200 independent runs, 0.0302. With naive error bars, only 17 % of runs would "contain" the true value instead of 68 %: the naive error bars claim 4.4 times more precision than the data have. Here, 20 000 samples are worth about 1000 independent ones.

Step 5

A real simulation: the pressure of the 9-atom box

Why: Now a quantity from an actual MD run: the instantaneous pressure of our 9-atom LJ box at T = 0.72, recorded every step for 100 000 steps. We compute the naive error, the block-averaging error, and check them against 20 independent runs.

What the numbers say: The honest result is P = 0.603 ± 0.036, ten times the naive ± 0.0035, and independent runs confirm it (0.0351). The system forgets its pressure only every 103 steps, so saving it every step gives no extra information beyond one sample per ~100 steps. The last rows show the limitation of block averaging: with 10 blocks the error estimate is itself uncertain by ±24 % and drifts up; read the plateau where there are still enough blocks. The uncertainty of an error from n blocks is about 1/√(2(n − 1)).

Step 6

Errors of derived quantities

Why: Results are often ratios or differences of measured numbers: a Haven ratio, a free-energy difference, a change of D with temperature. Independent errors add in squares: absolute errors for sums and differences, relative errors for products and ratios.

What the numbers say: The error of a ratio is dominated by its noisiest part: the 2.4 % of D_self hardly matters next to the 26 % of the collective term, so the Haven ratio inherits ±26 %. Differences are worse: two numbers known to 4–5 % give a difference known only to 28 %, which is why small free-energy differences or small changes between two force fields need long, careful runs. When the parts are correlated (computed from the same trajectory), these formulas are only approximate; the safest error bar then comes from repeating the whole analysis on independent runs or blocks.

Step 7

Equilibration and how long to run

Why: Two practical questions remain. If the run starts away from equilibrium (a box that is too large, velocities too hot), the start biases the average. And before running, one wants to know how many steps a target precision needs.

What the numbers say: Including the equilibration period shifts the mean by 0.30, twice the error bar, a systematic error that no amount of averaging removes. Discarding it costs a little precision (0.141 → 0.149) and removes the bias. Look at block means or running averages, and discard everything before they settle. For planning, the formula N_eff = (s/target)² × g turns a target precision into a number of steps: halving the error costs four times the run. These estimates assume the run samples all relevant states; slow events (Lab 4.6 dihedrals, Lab 7.4 barriers) need far more than the formula suggests.

Summary