MCMC Lecturenotes II Contents
1. Statistical Errors of Markov Chain MC Data 2. Autocorrelations 3. Integrated Autocorrelation Time and Binning 4. Illustration: Metropolis generation of normally distributed data 5. Self-consistent versus reasonable error analysis 6. Comparison of Markov chain MC algorithms
1
Statistical Errors of Markov Chain MC Data
In large scale MC simulation it may take months, possibly years, of computer time to collect the necessary statistics. For such data a thorough error analysis is a must. A typical MC simulation falls into two parts: 1. Equilibration: Initial sweeps are performed to reach the equilibrium distribution. During these sweeps measurements are either not taken at all or have to be discarded when calculating equilibrium expectation values. 2. Production: Sweeps with measurements are performed. Equilibrium expectation values are calculated from this statistics. A rule of thumb is: Do not spend more than 50% of your CPU time on measurements! The reason for this rule is that, that one cannot be off by a factor worse than two ( 2 in the statistical error). How many sweeps should be discarded for reaching equilibrium?
2
In a few situations this question can be rigorously answered with the Coupling from the Past method (Propp and Wilson). The next best thing to do is to measure the integrated autocorrelation time self-consistently and to discard, after reaching a visually satisfactory situations, a number of sweeps which is larger than the integrated autocorrelation time. In practice even this can often not be achieved. Therefore, it is re-assuring that it is sufficient to pick the number of discarded sweeps approximately right. With an increasing statistics the contribution of the non-equilibrium data dies out like 1/N , where N is the number of measurements. For large N effect is eventually swallowed by the statistical error, which declines the only like 1/ N . The point of discarding configurations for reaching equilibrium is that the factor in front of 1/N can be large. There can be far more involved situations, like that the Markov chain may end up in metastable configurations, which may even stay unnoticed (e.g. complex systems like spin glasses or proteins).
3
Autocorrelations
We like to estimate the expectation value f of some physical observable. We assume that the system has reached equilibrium. How many MC sweeps are needed to estimate f with some desired accuracy? To answer this question, one has to understand the autocorrelations within the Markov chain. Given is a time series of N measurements fi = fi(xi), i = 1, . . . , N from a Markov process, where xi are the configurations generated. The label i = 1, . . . , N runs in the temporal order of the Markov chain and the elapsed time, measured in updates or sweeps, between subsequent measurements fi, fi+1 is always the same, independently of i. The estimator of the expectation value f is f= 1 N fi .
4
With the notation t = |i - j| the definition of the autocorrelation function of the observable f is C(t) = Cij = (fi - fi ) (fj - fj ) = fifj - fi fj = f0ft - f 2 (1)
where we used that translation invariance in time holds for the equilibrium ensemble. The asymptotic behavior for large t is C(t) exp - t exp for t ,
where exp is called (exponential) autocorrelation time and is related to the second largest eigenvalue 1 of the transition matrix by exp = - ln 1
5
under the assumption that f has a non-zero projection on the corresponding eigenstate. Superselection rules are possible so that different autocorrelation times reign for different operators. The variance of f is a special case of the autocorrelations (1) C(0) = 2(f ) . Some algebra shows that the variance of the estimator f of the mean is now 2(f ) 2(f ) = 1+2 N
N -1
t=1
t 1- N
c(t)
with c(t) =
C(t) C(0)
.
This equation ought to be compared with the corresponding equation for uncorrelated random variables 2(f ) = 2(f )/N . The difference is the factor
6
in the bracket which defines the integrated autocorrelation time
N -1
int =
1+2
t=1
1-
t N
c(t)
.
For correlated data the variance of the mean is by the factor int larger than the corresponding naive variance for uncorrelated data: int 2(f ) 2(f ) 2 . = 2 with naive = N naive(f ) (2)
In most simulations one is interested in the limit N and int becomes
int = 1 + 2
t=1
c(t) .
(3)
The numerical estimation of the integrated autocorrelation time faces difficulties.
7
The variance of the estimator of int diverges for N :
int = 1 + 2
t=1
c(t) and 2( int)
because for large t each c(t) adds a constant amount of noise, whereas the signal dies out like exp(-t/exp). To obtain an estimate one considers the t-dependent estimator
t
int(t) = 1 + 2
t =1
c(t )
of the integrated autocorrelation time and looks out for a window in t for which int(t) is flat. To give a simple example, let us assume that the autocorrelation function is governed by a single exponential autocorrelation time C(t) = const exp - t exp .
8
In this case we can carry out the sum (3) for the integrated autocorrelation function and find
int = 1 + 2
t=1
e
-t/exp
2 e-1/exp . =1+ -1/exp 1-e
In particular, the difference between the asymptotic value and the finite t definition becomes then
int - int(t) = 2 e-t/exp
t =1
e-t /exp
2 e-(t+1)/exp . = -1/exp 1-e 1 the approximation
For a large exponential autocorrelation time exp int holds.
2 e-1/exp 2 - 2/exp =1+ = 2 exp - 1 2 exp =1+ = 1/exp 1 - e-1/exp
9
Integrated Autocorrelation Time and Binning
Using binning the integrated autocorrelation time can also be estimated via the variance ratio. We bin the time series into Nbs N bins of Nb = NBIN = NDAT N = Nbs NBINS
data each. Here [.] stands for Fortran integer division, i.e., Nb = NBIN is the largest integer N/Nbs, implying Nba Nb N . It is convenient to choose the values of N and Nbs so that N is a multiple of Nbs. The binned data are the averages jNb 1 N fi for j = 1, . . . , Nbs . fj b = Nb
i=1+(j-1)Nb
For Nb > exp the autocorrelations are essentially reduced to those between neighbor bins and even these approach zero under further increase of the binlength.
10
For a set of Nbs binned data fj b , (j = 1, . . . , Nbs) we may calculate the mean with its naive error bar. Assuming for the moment an infinite time series, we find the integrated autocorrelation time (2) from the following ratio of sample variances int =
Nb
b lim int
N
N
b with int =
N
s2 Nb f 2 sf
. (4)
In practice the Nb limit will be reached for a sufficiently large, finite value of Nb. The statistical error of the int estimate (4) is, in the first approximation, determined by the errors of s2 Nb . The typical situation is then that, due to the
f
central limit theorem, the binned data are approximately Gaussian, so that the error of s2 Nb is analytically known from the 2 distribution. Finally, the fluctuations of
f
s2 f
of the denominator give rise to a small correction which can be worked out.
Numerically most accurate estimates of int are obtained for the finite binlength Nb which is just large enough that the binned data are practically uncorrelated.
11
For applications it is convenient to choose N and Nb to be powers of 2. In the following we assume N = 2K , K 4 and Nb = 2Kb with Kb = 0, 1, . . . , K - 5, K - 4. (5)
Choosing the maximum value of Kb to be K - 4 implies that the smallest number of bins is min Nbs = 24 = 16 . (6) While the Student distribution shows that the confidence intervals of the error bars from 16 uncorrelated normal data are reasonable approximations to those of the Gaussian standard deviation, about 1000 independent data are needed to provide a decent estimate of the corresponding variance (at the 95% confidence level with an accuracy of slightly better than 10%). It makes sense to work with error bars from 16 binned data, but the error of the error bar, and hence a reliable estimate of int, requires far more data.
12
Illustration: Metropolis generation of normally distributed data
We generate normally distributed data according to the Markov process x = x + 2 a xr - a (7)
where x is the event at hand, xr a uniformly distributed random number in the range [0, 1), and the real number a > 0 is a parameter which relates to the efficiency of the algorithm. The new event x is accepted with the Metropolis probability Paccept(x ) = 1 for x 2 x2; exp[-(x 2 - x2)/2] for x
2
> x2.
(8)
If x is rejected, the event x is counted again. The Metropolis process an introduces autocorrelation time in the generation of normally distributed random data. We work with K = 17, i.e., N = 217 = 131072, data and take a = 3 for the Markov process (7), what gives an acceptance rate of approximately 50%.
13
The autocorrelation function (assignment a0401 01):
1 0.08 - 0.8 0.6 C(t) | 0.00 - 0.4 0.2 - 0.08 0 0 5 10 15 t 20 25 30
Figure 1: The autocorrelation function. Upper data: Metropolis time series for the normal distribution. Lower data: Gaussian random number generator. For t 11 the inlay shows the autocorrelations on an enlarged ordinate. The straight lines between the data points are just to guide the eyes. The curves start with C(0) 1 because the variance of the normal distribution is one.
14
Integrated autocorrelation time (assignment a0401 02):
4 3.5 3 int 2.5 2 1.5 1 0 20 1- 40 60 t 3- 2- 4- 3- 2- 1- 80 5 | 10 | 100 15 | 120 4- 10 | 20 | 30 |
Figure 2: Upper curves in the figure and its inlays: int(t). Lowest curve: Gaussian random number generator. Remaining curves: Binning estimators of the integrated autocorrelation time with one standard deviation bounds. The main figure relies on 221 data, the first inlay on 217 data, and the second inlay relies on 214 data.
15
b We compare the int estimators with the direct estimators int(t) at
N
t = Nb - 1 .
(9)
With this relation the estimators agree for binlength Nb = 1 and for larger Nb the relation gives the range over which we combine data into either one of the estimators. The approach of the binning procedure towards the asymptotic int value is slower than that of the direct estimate of int. For our large NDAT = 221 data set int(t) reaches its plateau before t = 20. All the error bars within the plateau are strongly correlated. Therefore, it is not recommendable to make an attempt to combine them. Instead, it is save to pick an appropriate single value and its error bar as the final estimate: int = int(20) = 3.962 0.024 from 221 = 2, 097, 152 data.
N
(10)
b The binning procedure, on the other hand, shows an increase of int all the way to
16
Nb = 27 = 128, where the estimate with the one confidence level error bounds is
128 3.85 int 3.94 from 214 = 16, 384 bins from 221 data .
How many data are needed to allow for a meaningful estimate of the integrated autocorrelation time? For the statistics of NDAT = 217 the autocorrelation signal disappears for t 11 into the statistical noise. Still, there is clear evidence of the hoped for window of almost constant estimates. A conservative choice is to take t = 20 again, which now gives int = int(20) = 3.86 0.11 from 217 data . Worse is the binning estimate, which for the 217 data is
32 3.55 int 3.71 from 212 = 4, 096 bins from 217 = 131, 072 data .
Our best value (10) is no longer covered by the two standard deviation zone.
17
For the second inlay the statistics is reduced to NDAT = 214. With the integrated autocorrelation time rounded to 4, this is 4096 times int. For binlength Nb = 24 = 16 we are then down to Nbs = 1024 bins, which are needed for accurate error bars of the error. To work with this number we limit, in accordance with equation (9), our int(t) plot to the range t 15. Still, we find a quite nice window of nearly constant int(t), namely all the way from t = 4 to t = 15. By a statistical fluctuation (assignment a0401 03) int(t) takes its maximum value at t = 7 and this makes int(7) = 3.54 0.13 a natural candidate. However, this value is inconsistent with our best estimate (10). The true int(t) increases monotonically as function of t, so we know that the estimators have become bad for t > 7. The error bar at t = 7 is obviously too small to take care of our difficulties. One may combine the t = 15 error bar. In this way the result is int = 3.54 0.21 for 214 = 16, 384 data, (11)
which achieves consistency in the two error bar range. For binlength Nb = 16 the
18
binning estimate is
16 2.93 int 3.20 from 210 = 1, 024 bins from 214 data.
(12)
Clearly, the binlength Nb = 16 is too small for an estimate of the integrated autocorrelation time. We learn from this investigation that one needs a binlength of at least ten times the integrated autocorrelation time int, whereas for the direct estimate it is sufficient to have t about four times larger than int.
19
Self-consistent versus reasonable error analysis
By visual inspection of the time series, one may get an impression about the length of the out-of-equilibrium part of the simulation. On top of this one should allow > int sweeps for the system to settle. A second reason why it appears necessary to control the the integrated autocorrelation times are the statistical errors of our measurements. Ideally the error bars are calculated as f = 2(f ) 2(f ) with (f ) = int . N
2
This constitutes a self-consistent error analysis of a MC simulation. However, the calculation of the integrated autocorrelation time may be out of reach. According to the Student distribution about twenty independent data are sufficient to estimate mean values with reasonably reliable error bars, while one thousand are needed for a 10%...