
% BB Oct 27 2006.
% BB Nov 21 2006.
% BB Feb 14 2004.

\documentclass[20pt,landscape,footrule]{foils}
\usepackage[dvips]{graphicx}
\usepackage[pdftex]{geometry}
\geometry{headsep=3ex,hscale=0.9}
\usepackage[pdftex]{color} % for the colored block
\date{\color{blue} Singapore MCMC Meeting, March 2004}
\usepackage[ignore]{pause}
%\usepackage{pause}
\usepackage{background}
%\usepackage{pp4slide}
%\MyLogo{AV HEP,CSIT}
\LogoOff
\begin{document}
\def\r{\color{red}}
\def\b{\color{blue}}
\def\s{\color{black}}

\foilhead{\r Lecturenotes 7 MCMC III -- Contents }\b \vskip 0.0pt
\begin{enumerate}
\item Statistical Errors of Markov Chain MC Data 
\item Autocorrelations
\item Integrated Autocorrelation Time and Binning
\item Illustration: Metropolis generation of normally
      distributed data
\item Self-consistent versus reasonable error analysis
\item Comparison of Markov chain MC algorithms
\end{enumerate}
 
\foilhead{\r Statistical Errors of Markov Chain MC Data } 
\vskip -30pt

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: 

\begin{enumerate}

{\r\item 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.

{\r\item Production:}\index{production: see data production}
\index{data production|textbf}Sweeps with measurements are performed.
Equilibrium expectation values are calculated from this statistics.

\end{enumerate}

A rule of thumb is: {\r 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 ($\sqrt{2}$ in the statistical error).

{\r How many sweeps should be discarded for reaching equilibrium?}

\newpage
In a few situations this question can be rigorously answered with 
the {\r 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 $\s 1/N$, 
where $\s N$ is the number of measurements. For large $\s N$ the effect 
is eventually swallowed by the statistical error, which declines only 
like $\s 1/\sqrt{N}$. The point of discarding configurations for 
reaching equilibrium is that the factor in front of $\s 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).
 
\foilhead{\r Autocorrelations} \vskip -30pt

We like to estimate the expectation value $\s \widehat{f}$ of some 
physical observable. We assume that the system has reached equilibrium.  
How many MC sweeps are needed to estimate $\s\widehat{f}$ with some 
desired accuracy? To answer this question, one has to understand the 
autocorrelations within the Markov chain.

Given is a {\r time series} of $\s N$ measurements 
$$\s f_i = f_i(x_i), ~~ i=1,\dots,N $$
from a Markov process, where $\s x_i$ are the configurations generated.
The label $\s i=1,\dots,N$ runs in the temporal order of the Markov 
chain and the elapsed time, measured in updates or sweeps, between 
subsequent measurements $\s f_i$, $\s f_{i+1}$ is always the same, 
independently of $\s i$.  The estimator of the expectation value 
$\s \widehat{f}$ is
$$\s \overline{f} = {1\over N} \sum f_i\ .$$
With the notation 
$$\s t=|i-j|$$ 
the definition of the {\r autocorrelation function} of the observable 
$\s \widehat{f}$ is 
\begin{equation} \label{autocorrelation} \s
 \widehat{C}(t) = \widehat{C}_{ij} = \langle\, (f_i - \langle f_i \rangle)\, 
 (f_j - \langle f_j\rangle)\, \rangle =
 \langle f_i f_j\rangle - \langle f_i \rangle\, \langle f_j\rangle = 
 \langle f_0 f_t\rangle - \widehat{f}^{\,2} \end{equation}
where we used that translation invariance in time holds for the 
equilibrium ensemble. The asymptotic behavior for large $t$ is
$$\s \widehat{C}(t) \sim \exp \left( - {t\over \tau_{\rm exp}} \right)
 ~~{\rm for}~~ t \to \infty , $$
where $\s \tau_{\rm exp}$ is called {\r (exponential) autocorrelation 
time} and is related to the second largest eigenvalue $\s \lambda_1$ of 
the transition matrix by
$$\s \tau_{\rm exp} = - \frac{1}{\ln\lambda_1} $$
under the assumption that $\s 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 $\s f$ is a special case of the 
autocorrelations~(\ref{autocorrelation})
$$\s \widehat{C} (0) = \sigma^2 (f)\ .$$
Some algebra shows that the variance of the estimator 
$\s \overline{f}$ of the mean is now 
$$\s \sigma^2(\overline{f})\ =\ {\sigma^2(f)\over N} \left[ 1 +
  2 \sum_{t=1}^{N-1}\left( 1-{t\over N}\right)\,\widehat{c}(t)\right]~~
  {\rm with}~~\widehat{c}(t) = \frac{\widehat{C}(t)}{\widehat{C}(0)}\ .$$
This equation ought to be compared with the corresponding equation 
for uncorrelated random variables $\s \sigma^2(\overline{f}) = 
{\sigma^2(f)/ N}$. The difference is the factor in the bracket which
defines the {\r integrated autocorrelation time}
$$\s \tau_{\rm int}\ =\ \left[ 1 + 2 \sum_{t=1}^{N-1} 
 \left(1-{t\over N}\right)\, \widehat{c}(t) \right]\ .$$
For correlated data the variance of the mean is by the factor 
$\s \tau_{\rm int}$ larger than the corresponding
{\r naive variance} for uncorrelated data:
\begin{equation} \label{variance_ratio} \s
  \tau_{\rm int}\ =\ {\sigma^2(\overline{f})\over 
  \sigma^2_{\rm naive} (\overline{f}) } ~~{\rm with}~~
  \sigma^2_{\rm naive} = {\sigma^2(f)\over N}\ .
\end{equation}
In most simulations one is interested in the limit $N\to\infty$ and 
$\s\tau_{\rm int}$ becomes
\begin{equation} \label{ia_time} \s
  \tau_{\rm int}\ =\ 1 + 2 \sum_{t=1}^{\infty} \widehat{c}(t)\ . 
\end{equation}
The numerical estimation of the integrated autocorrelation time faces 
difficulties. The variance of the estimator of $\s\tau_{\rm int}$ 
diverges for $\s N\to\infty$:
$$\s \overline{\tau}_{\rm int}\ =\ 1 + 2 \sum_{t=1}^{\infty}
     \overline{c}(t) ~~{\rm and}~~
     \sigma^2 (\overline{\tau}_{\rm int})\ \to\ \infty $$
because for large $\s t$ each $\s \overline{c}(t)$ adds a constant 
amount of noise, whereas the signal dies out like 
$\s \exp (-t/\tau_{\rm exp})$. To obtain an estimate one considers 
the $\s t$-dependent estimator
$$\s \overline{\tau}_{\rm int}(t)\ =\ 1 + 2 \sum_{t'=1}^t 
     \overline{c}(t') $$
of the integrated autocorrelation time and looks out for a {\r window} 
in $\s t$ for which $\s \overline{\tau}_{\rm int}(t)$ is flat. To give 
a {\bf simple example}, let us assume that the autocorrelation function 
is governed by a single exponential autocorrelation time
$$\s\widehat{C}(t)\ =\ 
    const\,\exp\left(-{t\over\tau_{\rm exp}}\right)\ .$$
In this case we can carry out the sum~(\ref{ia_time}) for the 
integrated autocorrelation function and find
$$\s \tau_{\rm int} = 1 + 2 \sum_{t=1}^{\infty}
     e^{-t/\tau_{\rm exp}} = 1 + { 2\,e^{-1/ \tau_{\rm exp}} 
     \over 1 - e^{-1/ \tau_{\rm exp} } }\ .$$
In particular, the difference between the asymptotic value and the 
finite $\s t$ definition becomes then
$$\s \tau_{\rm int} - \tau_{\rm int}(t) = 2\, e^{-t/\tau_{\rm exp}} 
     \sum_{t'=1}^{\infty} e^{-t'/\tau_{\rm exp}} = 
     {2\,e^{-(t+1)/\tau_{\rm exp}}\over 1-e^{-1/\tau_{\rm exp}} }\ .$$
For a large exponential autocorrelation time $\s \tau_{\rm exp}\gg 1$
the approximation
$$\s \tau_{\rm int} = 1 + { 2\,e^{-1/ \tau_{\rm exp}} \over
     1 - e^{-1/ \tau_{\rm exp} } } \cong
     1 + { 2 - 2 / \tau_{\rm exp} \over 1 / \tau_{\rm exp} } =
     2\,\tau_{\rm exp} - 1\ \cong\ 2\,\tau_{\rm exp} $$
holds. 

\foilhead{\r Integrated Autocorrelation Time and Binning}\vskip -30pt

Using binning the integrated autocorrelation time can also be estimated 
via the variance ratio.  We bin the time series into $\s N_{bs}\le N$ 
bins of
$$\s N_b = {\tt NBIN} = \left[ N\over N_{bs} \right] =
     \left[ {\tt NDAT} \over {\tt NBINS} \right] $$
data each. Here $\s [.]$ stands for Fortran integer division, {\it 
i.e.}, $\s N_b={\tt NBIN}$ is the largest integer $\s \le N/N_{bs}$, 
implying $\s N_{ba}\cdot N_b\le N$. It is convenient to choose the 
values of $\s N$ and $\s N_{bs}$ so that $\s N$ is a multiple of 
$\s N_{bs}$. The binned data are the averages
$$\s f^{N_b}_j = {1\over N_b} \sum_{i=1+(j-1)N_b}^{jN_b} f_i 
     ~~~~{\rm for}~~ j=1,\dots ,N_{bs}\ .$$
For $\s N_b > \tau_{\rm exp}$ the autocorrelations are essentially 
reduced to those between neighbor bins and even these approach zero 
under further increase of the binlength.

% \newpage
For a set of $\s N_{bs}$ binned data $\s f^{N_b}_j$, $\s (j=1,\dots ,
N_{bs})$ we may calculate the mean with its {\r naive error bar}. 
Assuming for the moment an infinite time series, we find the integrated 
autocorrelation time~(\ref{variance_ratio}) from the following ratio of 
sample variances
\begin{equation} \label{iat_binning} \s
\tau_{\rm int}\ =\ \lim_{N_b\to\infty}\, \tau_{\rm int}^{N_b}\ 
~~{\rm with}~~\ \tau_{\rm int}^{N_b}\ =\ \left( 
{ s^2_{\overline{f}^{N_b}} \over s^2_{\overline{f}} }\right)\ .
\end{equation}
In practice the $\s N_b\to\infty$ limit will be reached for a 
sufficiently large, finite value of $\s N_b$. The statistical error 
of the $\s \tau_{\rm int}$ estimate (\ref{iat_binning}) is, in the 
first approximation, determined by the errors of $\s 
s^2_{\overline{f}^{N_b}}$.  The typical situation is then that, 
due to the central limit theorem, the binned data are approximately 
Gaussian, so that the {\r error of $\s s^2_{\overline{f}^{N_b}}$ is 
analytically known} from the $\s \chi^2$ distribution. Finally, the 
fluctuations of $\s s^2_{\overline{f}}$ of the denominator give 
rise to a small correction which can be worked out. 

Numerically most accurate estimates of $\s \tau_{\rm int}$ are obtained 
for the finite binlength $\s N_b$ which is just large enough that the 
binned data are practically uncorrelated. 

\newpage
For applications it is convenient to choose $\s N$ and $\s N_{b}$ to be 
powers of {\s 2}. In the following we assume
\begin{equation}\s \label{num_data}
  N = 2^K,\ K\ge 4 ~~{\rm and}~~ N_b= 2^{K_b} ~~{\rm with}~~
  K_b=0,1,\dots ,K-5,K-4 .
\end{equation}
Choosing the maximum value of $\s K_b$ to be $\s K-4$ implies that the 
smallest number of bins is 
\begin{equation}\s \label{iat_error}
  N_{bs}^{\min}\ =\ 2^4\ =\ 16\ .
\end{equation}
While the Student distribution shows that the confidence intervals 
of the error bars from {\s 16} uncorrelated normal data are reasonable 
approximations to those of the Gaussian standard deviation,
about {\s 1000} independent data are needed to provide a decent 
estimate of the corresponding variance (at the {\s 95\%} confidence 
level with an accuracy of slightly better than {\s 10\%}). 

It makes sense to work with error bars from {\s 16} binned data, but 
the error of the error bar, and hence a reliable estimate of 
$\s \tau_{\rm int}$, requires far more data.
 
\foilhead{\r Illustration: Metropolis generation of normally
             distributed data} \vskip -30pt

We generate normally distributed data according to the Markov process
\begin{equation}\s \label{gau_metro}
  x' = x + 2\,a\,x^r - a
\end{equation}
where $\s x$ is the event at hand, $\s x^r$ a uniformly distributed 
random number in the range $\s [0,1)$, and the real number $\s a>0$ 
is a parameter which relates to the efficiency of the algorithm. The 
new event $\s x'$ is accepted with the Metropolis probability
\begin{equation}\s \label{mg_accept}
 P_{\rm accept}(x') = \cases{ 1~~{\rm for}~~x'^{\,2}\le x^2;\cr
  \exp [-( x'^{\,2}-x^2)/2 ]\ {\rm for}\ x'^{\,2}>x^2 .}
\end{equation}
If $\s x'$ is rejected, the event $\s x$ is counted again. The 
Metropolis process introduces an autocorrelation time in the 
generation of normally distributed random data.

We work with $\s K=17$, {\it i.e.}, $\s N=2^{17}=131072$, data and take
$\s a=3$ for the Markov process (\ref{gau_metro}), what gives 
an acceptance rate of approximately {\s 50\%}. 

\newpage
{\r The autocorrelation function (assignment {\tt a0401\_01}):}

\begin{figure} \begin{center}
\includegraphics[height=0.59\textheight]{figures/acor.pdf}
\caption{The autocorrelation function. Upper data: Metropolis time 
series for the normal distribution. Lower data: Gaussian random number 
generator. For $t\ge 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 $\overline{C}(0)\approx 1$ 
because the variance of the normal distribution is one.}
\end{center} \end{figure}

\newpage
{\r Integrated autocorrelation time (assignment {\tt a0401\_02}):}

\begin{figure} \begin{center}
\includegraphics[height=0.60\textheight]{figures/at_int.pdf}
\caption{Upper curves in the figure and its inlays: $\tau_{\rm 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 $2^{21}$ data, 
the first inlay on $2^{17}$ data, and the second inlay relies on 
$2^{14}$ data.}
\end{center} \end{figure}

\newpage
We compare the $\s \tau^{N_b}_{\rm int}$ estimators with the direct 
estimators $\s \tau_{\rm int}(t)$ at
\begin{equation}\s \label{t_eq_Nbm1}
  t = N_b -1\ .
\end{equation}
With this relation the estimators agree for binlength $\s N_b=1$ and for 
larger $\s N_b$ 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 $\s \tau_{\rm int}$ value is slower than that of 
the direct estimate of $\s \tau_{\rm int}$.

For our large $\s {\tt NDAT}=2^{21}$ data set $\s \tau_{int}(t)$ reaches 
its plateau before $\s 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:
\begin{equation}\s \label{iat_Gauss_estimate}
  \tau_{\rm int} = \tau_{\rm int} (20) = 3.962 \pm 0.024 ~~{\rm from}~~
  2^{21}=2,097,152\ {\rm data}.
\end{equation}
The binning procedure, on the other hand, shows an increase of 
$\s \tau^{N_b}_{\rm int}$ all the way to $\s N_b=2^7=128$, where 
the estimate with the one confidence level error bounds is
$$\s  3.85 \le \tau^{128}_{\rm int} \le 3.94 ~~{\rm from}~~
      2^{14}=16,384\ {\rm bins}~~{\rm from}~~2^{21}\ {\rm data}\ .$$
% \newpage
{\r How many data are needed to allow 
for a meaningful estimate of the integrated autocorrelation time?}

For the statistics of $\s {\tt NDAT}=2^{17}$ the autocorrelation signal 
disappears for $\s t\ge 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 $\s t=20$ again, which now gives
$$\s  \tau_{\rm int} = \tau_{\rm int}(20)=3.86\pm 0.11
    ~~{\rm from}~~2^{17}\ {\rm data}\ . $$
Worse is the binning estimate, which for the $\s 2^{17}$  data is 
$$\s 3.55 \le \tau^{32}_{\rm int} \le 3.71 ~~{\rm from}~~
     2^{12}=4,096\ {\rm bins}~~{\rm from}~~2^{17}=131,072\ 
     {\rm data}\ . $$
Our best value (\ref{iat_Gauss_estimate}) is no longer covered by 
the two standard deviation zone.

For the second inlay the statistics is reduced to $\s {\tt NDAT}=
2^{14}$. With the integrated autocorrelation time rounded to {\s 4}, 
this is {\s 4096} times $\s \tau_{\rm int}$. For binlength $\s N_b=2^4
=16$ we are then down to $\s N_{bs}=1024$ bins, which are needed for 
accurate error bars of the error. To work with this number we limit, 
in accordance with equation (\ref{t_eq_Nbm1}), our $\s 
\tau_{\rm int}(t)$ plot to the range $\s t\le 15$. Still, we find a 
quite nice window of nearly constant $\s \tau_{\rm int}(t)$, namely 
all the way from $\s t=4$ to $\s t=15$. By a statistical fluctuation
(assignment {\tt a0401\_03}) $\s \tau_{\rm int}(t)$ takes its maximum 
value at $\s t=7$ and this makes $\s \tau_{\rm int}(7) = 3.54\pm 0.13$ 
a natural candidate. However, this value is inconsistent with our best 
estimate (\ref{iat_Gauss_estimate}). The true $\s \tau_{\rm int}(t)$ 
increases monotonically as function of $\s t$, so we know that the 
estimators have become bad for $\s t>7$. The error bar at $\s t=7$ 
is obviously too small to take care of our difficulties. One may 
combine the $\s t=15$ error bar. In this way the result is 
\begin{equation}\s  \label{tau_int14}
  \tau_{\rm int}=3.54\pm 0.21~~{\rm for}~~2^{14}=16,384\ {\rm data},
\end{equation}
which achieves consistency in the two error bar range. For binlength 
$\s N_b=16$ the binning estimate is 
\begin{equation}\s \label{F_cl_example}
  2.93 \le \tau^{16}_{\rm int} \le 3.20 ~~{\rm from}~~
  2^{10}=1,024\ {\rm bins}~~{\rm from}~~2^{14}\ {\rm data}. 
\end{equation}
Clearly, the binlength $\s N_b=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 $\s \tau_{\rm int}$, whereas for the direct 
estimate it is sufficient to have $\s t$ about four times larger 
than $\s \tau_{\rm int}$.
 
\foilhead{\r Self-consistent versus reasonable error analysis}
\vskip -30pt
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 $\s >\tau_{\rm 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
$$\s \triangle \overline{f}\ =\ \sqrt{\sigma^2(\overline{f})}
     ~~~{\rm with}~~~ \sigma^2(\overline{f})\ =\ 
     \tau_{\rm int}\,{\sigma^2(f)\over N}\ .$$
This constitutes a {\r 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 {\r twenty}
independent data are sufficient to estimate mean values with reasonably 
reliable error bars, while {\r one thousand} are needed for a {\s 10\%}
error on the error bar at the two sigma confidence level.

{\r In practice, one has to be content with what can be done.} Often 
this means to rely on the binning method. We simply calculate error 
bars of our ever increasing statistics with respect to a fixed number 
of {\s \tt NBINS} bins, where 
$$\s {\tt NBINS}\ \ge\ 16$$
are large enough values. In addition, we may put {\s 10\%} of the 
initially planned simulation time apart for reaching equilibrium. 
{\it A-posteriori}, this can always be increased. Once the statistics 
is large enough, this means at least $\s t_{\max} > {\tt NBINS}\, 
\tau_{\rm int}$ for the largest $\s t$ of our time series, our small 
number of binned data become effectively independent.

How do we know that the statistics has become large enough? In 
practical applications there can be indirect arguments, which tell 
or suggest us that the integrated autocorrelation time is in fact 
(much) smaller than the achieved bin length. This is no longer 
self-consistent, as we perform no explicit measurement of 
$\s \tau_{\rm int}$, but it is a {\r reasonable error analysis}.
 
\foilhead{\r Comparison of Markov chain MC algorithms}
\vskip -30pt
Is the 1-hit Metropolis algorithm more efficient with sequential 
updating or with random updating (assignment {\tt a0402\_01} B).
%
\begin{figure} \begin{center}
\includegraphics[height=0.52\textheight]{figures/iat1r.pdf}
\caption{Integrated autocorrelation times of the Metropolis process 
with random updating versus sequential updating for the $d=2$ Ising 
model at $\beta =0.4$ The ordering of the curves is identical with 
that of the labels in the figure.}
\end{center} \end{figure}

\newpage
{\r The $\s d=2$ Ising model off the critical point  
(assignment {\tt a0402\_01}~A).}
% {\r The $\s d=2$ Ising model off and on the critical point
% (assignment {\tt a0402\_01}~A).}
%
\begin{figure} \begin{center}
\includegraphics[height=0.60\textheight]{figures/iat1.pdf}
\caption{One-hit Metropolis algorithm with sequential updating: Lattice 
size dependence of the integrated autocorrelation time for the $d=2$ 
Ising model at $\beta =0.4$. \r The ordering of the curves is identical 
with that of the labels in the figure.}
\end{center} \end{figure}

\newpage
{\r Another MC dynamics, Swendsen-Wang (SW) and Wolff (W) cluster 
algorithm:}

\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/at_cl.pdf}
\caption{Estimates of integrated autocorrelation times from 
simulations of the $d=2$ Ising model at the critical temperature 
$\beta_c=0.44068679351$ (assignment {\tt a0503\_05}). }
\end{center} \end{figure}

\newpage
{\r Metropolis versus heat bath updating for the {\s 10-state
$d=2$} Potts model:}

\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/iathb.pdf}
\caption{Systematic updating: Comparison of the integrated 
autocorrelation times of the 1-hit and 2-hit Metropolis algorithms 
and the heat bath algorithm for the 10-state Potts model on $L\times L$ 
lattices at $\beta =0.62$ (assignment {\tt a0402\_06}). The $L=40$ 
and $L=80$ curves lie almost on top of one another. }
\end{center} \end{figure}

\end{document}

\newpage
{\r Critical slowing down:} Increase of the $\s \tau_{\rm int}\sim L^z$ 
at the critical point, where $\s z\approx 2$ is the {\r dynamical 
critical exponent} (assignment {\tt a0402\_02}~D). Estimates of $\s z$ 
are compiled in the book by Landau and Binder.
\begin{figure} \begin{center}
\includegraphics[height=0.54\textheight]{figures/iat2_large.pdf}
\caption{One-hit Metropolis algorithm with sequential updating: 
$L$ dependence of the integrated autocorrelation time for the $d=2$ 
Ising model at its critical temperature. The ordering of the 
curves is identical with that of the labels.}
\end{center} \end{figure}
