
% BB Oct 10 2008.
% BB Oct 01 2007.
% BB Nov 21 2006.
% BB Mar 11 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 4 MCMC I -- Contents}\b 
\begin{enumerate}
\item Statistical Physics and Potts Models
\item Sampling and Re-weighting
\item Importance Sampling and Markov Chain Monte Carlo
\item The Metropolis Algorithm 
\item The Heatbath Algorithm (Gibbs Sampler)
\item Start and Equilibration
\item Energy Checks
\item Specific Heat
\end{enumerate} \vskip -30pt
 
\foilhead{\r Statistical Physics and Potts Model} \vskip -30pt

MC simulations of systems described by the Gibbs canonical ensemble 
aim at calculating estimators of physical observables at a {\r 
temperature} {\s $T$.} In the following we choose units so that $\s
\beta = 1/T$ and consider the calculation of the {\r expectation 
value} of an {\r observable} $\s\mathcal{O}$. Mathematically all 
systems on a computer are discrete, because a finite word length 
has to be used. Hence, % the expectation value is given by the sum
{\s \begin{eqnarray} \label{O}
 \widehat{\mathcal{O}} = \widehat{\mathcal{O}} (\beta) = 
 \langle \mathcal{O} \rangle &=& Z^{-1} \sum_{k=1}^K
 \mathcal{O}^{(k)}\,e^{-\beta\,E^{(k)} }\\
{\rm where}~~~
 Z\ =\ Z(\beta) &=& \sum_{k=1}^K e^{-\beta\,E^{(k)} }
\end{eqnarray} }
is the {\r partition function}. The index $\s k=1,\dots , K$ labels 
all {\r configurations (or microstates)} of the system, and $\s
E^{(k)}$ is the (internal) {energy } of configuration $\s k$. To 
distinguish the configuration index from other indices, it is put 
in parenthesis.

\newpage
We introduce {\r generalized Potts models} in an  external magnetic 
field on $\s d$-dimensional hypercubic lattices with periodic boundary
conditions. Without being overly complicated, these models are general
enough to illustrate the essential features we are interested in.
In addition, various subcases of these models are by themselves of 
physical interest. Generalizations of the algorithmic concepts to 
other models are straightforward, although technical complications 
may arise. 

We define the energy function of the system by
{\s \begin{equation} \label{E}
 -\beta\,E^{(k)} = -\beta\,E_0^{(k)} + H\,M^{(k)} 
\end{equation} }
where
{\s \begin{equation} \label{E0}
  E_0^{(k)} = -2\,\sum_{\langle ij\rangle} J_{ij}(q_i^{(k)},q_j^{(k)})\
  \delta (q_i^{(k)},q_j^{(k)}) + {2\,d\,N\over q}
\end{equation}
$$ {\rm with}~~\ \delta (q_i,q_j)
 = \left\{ \begin{array}{c} 1\ {\rm for}\ q_i=q_j \\
         0\ {\rm for}\ q_i\ne q_j \end{array}  \right. 
~~~{\rm and}~~~
 M^{(k)} = 2 \sum_{i=1}^N \delta (1,q_i^{(k)})\, .  $$ }

\newpage
The sum $\s\langle ij\rangle$ is over the nearest neighbor lattice 
sites and $\s q_i^{(k)}$ is called the {\r Potts spin} or {\r Potts 
state} of configuration $\s k$ at site $\s i$. For the $\s q$-state 
Potts model $\s q^{(k)}_i$ takes on the values $\s 1,\dots ,q$. The 
external magnetic field is chosen to interact with the state $\s q_i=1$ 
at each site $i$, but not with the other states $\s q_i\ne 1$. The 
{\s $J_{ij} (q_i,q_j)$, $(q_i=1,\dots ,q;\, q_j=1,\dots, q)$} functions 
define the exchange coupling constants between the states at site $\s i$ 
and site~$\s j$.  The energy function describes a number of physically
interesting situations. With
\begin{equation} \label{Jij} \s
 J_{ij} (q_i,q_j) \equiv J > 0~~\ ({\rm conventionally}\ J=1) 
\end{equation}
the original model is recovered and $\s q=2$ becomes equivalent to the 
Ising ferromagnet. The Ising case of Edwards-Anderson spin glasses and 
quadrupolar Potts glasses are obtained when the exchange constants are 
quenched random variables. Other choices of the $\s J_{ij}$ include 
anti-ferromagnets and the fully frustrated Ising model.

For the {\r energy per spin} the notation is:\ $\s~e_s = E/N\,$.

\newpage
The normalization is chosen so that $e_s$ agrees for $\s q=2$ with
the conventional Ising model definition, $\s\beta = \beta^{\rm Ising} 
= \beta^{\rm Potts}/2\,$.

For the $\s 2d$ Potts models a number of exact results are known in 
the infinite volume limit, mainly due to work by Baxter. The phase 
transition temperatures are
\begin{equation} \label{2d_Potts_bc} \s
  {1\over 2}\, \beta_c^{\rm Potts}\ =\ \beta_c\ =\ {1\over T_c}\ 
  =\ {1\over 2}\, \ln (1 + \sqrt{q} ), ~~q=2,3,\dots\ .
\end{equation}
At $\s\beta_c$ the average energy per state is 
\begin{equation} \label{Potts_e0s} \s
 e_{0s}^c = E_0^c/N = {4\over q} - 2 - 2/\sqrt{q}\ .
% = {4\over q} + e_{0s}[{\rm Baxter}]
\end{equation}
The phase transition is second order for $\s q\le 4$ and first order 
for $\s q\ge 5$ for which the exact infinite volume {\r latent heats} 
$\s \triangle e_{0s}$ and {\r entropy jumps} $\s \triangle s$ were 
also found by Baxter, while the {interface tensions $\s f_s$ were 
derived later.

\foilhead{\r Some Potts Configurations \s ($2d$, $q=10$)}

\begin{enumerate}
\item Ordered with small fluctuations.
\item Disordered droplet in ordered background.
\item Percolated disordered phase (order-disorder separation)
\end{enumerate}

\newpage

\begin{figure} \begin{center}
\includegraphics[height=1.2\textheight]{figures/color_potts1.jpg}
% \caption{\small .  }
\end{center} \end{figure}

\begin{figure} \begin{center}
\includegraphics[height=1.2\textheight]{figures/color_potts2.jpg}
% \caption{\small .  }
\end{center} \end{figure}

\begin{figure} \begin{center}
\includegraphics[height=1.2\textheight]{figures/color_potts3.jpg}
% \caption{\small .  }
\end{center} \end{figure}

\foilhead{\r Sampling and Re-weighting}\vskip -30pt

For the Ising model it is straightforward to {\r sample statistically 
independent configurations}. We simply have to generate $\s N$ spins, 
each either up or down, each with 50\% likelihood. This is called {\r 
random sampling}, see the figure.
%
\begin{figure} \begin{center}
\includegraphics[height=0.4\textheight]{figures/2di_h_es.pdf}
\caption{\small 
Energy histograms of 100$\,$000 entries each for the Ising model on 
an $20\times 20$ lattice: Random Sampling gives statistically independent
configurations at $\beta=0$. Histograms at $\beta=0.2$ and $\beta=0.4$ 
are generated by Markov chain MC. Re-weighting of the $\beta=0$ random 
configurations to $\beta=0.2$ is shown to fail (assignments
{\tt a0301\_02} and {\tt a0303\_02}).  }
\end{center} \end{figure}
\newpage

Is is very important to {\r distinguish} the energy {\r measurements} 
$\s e_s$ on single configurations from the {\r expectation value}
expectation value $\s\widehat{e}_s$. The expectation value 
$\s\widehat{e}_s$ is a single number, while $\s e_s$ fluctuates. 
From the measurement of many $\s e_{s}$ values one finds {\r 
estimators} of its moments. The mean is is denoted by 
$\s\overline{e}_{s}$ and fluctuates.

Histogram entries at $\s \beta=0$ can be re-weighted, so that they 
correspond to other $\s \beta$ values. We simply have to multiply the 
entry corresponding to energy $\s E$ by $\s c_{\beta}\exp (-\beta E)$.
Similarly histograms corresponding to the Gibbs ensemble at some 
value $\s\beta_0$ can be re-weighted to other $\s\beta$ values. 
(Care has to be taken to ensure that the involved arguments of the 
exponential function do not become too large! This can be done by
logarithmic coding.) 

{\r Re-weighting} has a long history. For Finite Size Scaling (FSS) 
investigations of second order phase transitions its usefulness 
has been stressed by {\r Ferrenberg and Swendsen} (accurate 
determinations of peaks of the specific heat or of susceptibilities). 
\newpage

In the figure re-weighting is done from $\s\beta_0=0$ to
$\s\beta=0.2$. But, by comparison to the histogram from a Metropolis 
MC calculation at $\s\beta=0.2$, the result is seen to be disastrous. 
The reason is easily identified: In the range where the $\s\beta=0.2$ 
histogram takes on its maximum, the $\s\beta=0$ histogram  has not a 
single entry, {\it i.e.}, our naive sampling procedure misses the 
{\r important configurations} at $\s\beta =0.2$. Re-weighting to new 
$\s\beta$ values works only in a range $\s \beta_0\pm \triangle \beta$, 
where $\s\triangle\beta\to 0$ in the infinite volume limit. 

{\bf Important Configurations}

Let us determine the important contributions to the partition function. 
The partition function can be re-written as a sum over energies
\begin{equation} \label{Z_sum_E} \s
  Z = Z(\beta) = \sum_E n(E)\, e^{-\beta\,E}
\end{equation}
where the unnormalized spectral density $\s n(E)$ is defined as the 
number of microstates $\s k$ with energy $\s E$ (remember, on the 
computer energy values are always discrete). 
 
\newpage

For a fixed value of $\s \beta$ the energy probability density
\begin{equation} \label{P_E} \s
 P_{\beta} (E) = c_{\beta}\, n(E)\, e^{-\beta E} 
\end{equation}
is peaked around the average value $\s\widehat{E}(\beta)$, where 
$\s c_{\beta}$ is a normalization constant so that the $\s\sum_E 
P_{\beta}(E)=1$ holds (consider the limits $\s\beta=0$ and 
$\s\beta\to\infty.$)

Away from phase transitions the width of the 
energy distribution is $\s \triangle E\sim\sqrt{V}$. This follows from 
the fact that, away from phase transition points, the fluctuations of 
the $\s N\sim V$ lattice spins are essentially uncorrelated, so that 
the magnitude of a typical fluctuations is $\s\sim\sqrt{N}$. From this 
we find that the re-weighting range is $\s\triangle\beta\sim 1/\sqrt{V}$, 
as the energy is an extensive quantity $\s\sim V$ so that $\s\triangle 
\beta E\sim \sqrt{V}$ can stay within the fluctuation of the system. 

Interestingly, the re-weighting range increases at \r second order phase 
transitions, \b because \r critical fluctuations \b are larger than 
non-critical fluctuations. Namely, one has $\s \triangle E \sim V^x$ 
with $\s 1/2<x<1$ and the requirement $\s \triangle\beta E\sim V^x$ 
yields $\s \triangle\beta \sim V^{x-1}$. 
 
\newpage

For \r first order phase transitions \b one has a latent heat 
$\s\triangle V\sim V$, but this does not mean that the re-weighting 
range becomes of order one. In essence, the fluctuations collapse, 
because the two phases become separated by an interfacial tension. 
One is back to 
fluctuations within either of the two phases, {\it i.e.} 
$\s\triangle\beta\sim 1/\sqrt{V}$.

The important configurations at temperature $\s T=1/\beta$ are at the 
energy values for which the probability density $\s P_{\beta}(E)$ is 
large. To sample them efficiently, one needs a procedure which 
generates the configurations with their Boltzmann weights
\begin{equation} \label{w_B} \s
  w_B^{(k)} = e^{-\beta E^{(k)}}\ . 
\end{equation}
The number of configurations $\s n(E)$ and the weights combine then 
so that the probability to generate a configuration at energy $\s E$
becomes precisely $\s P_{\beta}(E)$ as given by equation~(\ref{P_E}).

\foilhead{\r Importance Sampling and Markov Chain Monte Carlo}
\vskip -30pt

For the canonical ensemble {\r importance sampling} generates 
configurations $\s k$
with probability
\begin{equation} \label{P_B} \s
  P_B^{(k)} = c_B\, w^{(k)}_B = c_B\, e^{-\beta E^{(k)}}\,,\qquad
  c_B~~{\rm constant}\,.
\end{equation}
The {\r state vector} $\s (P_B^{(k)})$, for which the configurations 
are the vector indices, is called {\r Boltzmann state}. 
The {\r expectation value} becomes the {\r arithmetic average:}
\begin{equation} \label{O_B} \s
  {\widehat {\cal O}} = {\widehat {\cal O}} (\beta) =
  \langle {\cal O} \rangle  = \lim_{N_K\to \infty} 
  {1\over N_K} \sum_{n=1}^{N_K} {\cal O}^{(k_n)}\ .
\end{equation}
When the sum is truncated we obtain an {\r estimator of the 
expectation value:}
\begin{equation} \label{estimate_O_B} \s
 {\overline {\cal O}} =
 {1\over N_K} \sum_{k=1}^{N_K} {\cal O}^{(k_n)}\ .
\end{equation}
% \newpage

Normally, we cannot generate configurations $\s k$ directly with
probability~(\ref{P_B}). But they may be found as members of the 
equilibrium distribution of a dynamic process. In practice {\r Markov 
chains} are used. A {\r Markov process} is a particularly simple dynamic 
process, which generates configuration $\s k_{n+1}$ stochastically from 
configuration $\s k_n,$ so that no information about previous 
configurations $\s k_{n-1}, k_{n-2},\,\dots$ is needed. The elements 
of the Markov chain {\r time series} are the configurations. Assume 
that the configuration $\s k$ is given. Let the {\r transition 
probability} to create the configuration $\s l$ in one step from 
$\s k$ be given by $\s W^{(l)(k)} = W[k\to l]$. In essence, the matrix
\begin{equation} \label{W} \s
  W = \left( W^{(l)(k)}\right)
\end{equation}
defines the Markov process. Note, that this {\r transition matrix}
is a very big (never stored in the computer!), because its labels 
are the configurations. To achieve our goal to generate configurations 
with the desired probabilities, the matrix $\s W$ is required 
to satisfy the following properties:
\newpage

\begin{description}

\item{(i)} \index{ergodicity}{\r\bf Ergodicity (Irreducibility in 
Math Literature):}
\begin{equation} \label{ergodicity} \s
  e^{-\beta E^{(k)}} > 0 ~~{\rm and}~~ e^{-\beta E^{(l)}} > 0
  ~~~{\rm imply:}
\end{equation}
an integer number $\s n>0$ exists so that $\s (W^n)^{(l)(k)}>0$ holds.

\item{(ii)} {\r\bf Normalization:}
\begin{equation} \label{normalization} \s
 \sum_l W^{(l)(k)} = 1\ .
\end{equation}

\item{(iii)} \index{balance|textbf}{\r\bf Balance (Stationarity):} 
\begin{equation} \label{balance} \s
 \sum_k W^{(l)(k)}\, e^{-\beta E^{(k)}}\ =\ e^{-\beta E^{(l)}}\ .
\end{equation}
Balance means: The Boltzmann state~(\ref{P_B}) is an eigenvector with 
eigenvalue $\s 1$ of the transition matrix $\s W = (W^{(l)(k)})$.

\end{description}
\newpage

In statistical physics the ensemble notation is frequently used. By 
definition, an {\r ensemble} is a collection of configurations so that 
to each configuration $\s k$ a probability $\s P^{(k)}$ is assigned, 
$\s \sum_k P^{(k)} = 1$. The {\r Gibbs or Boltzmann ensemble} $\s E_B$ 
is defined to be the ensemble with probability distribution~(\ref{P_B}). 

An {\r equilibrium ensemble} $\s E_{eq}$ of the Markov process is 
defined by its probability distribution $\s P_{eq}$ satisfying
\begin{equation} \label{E_eq} \s
 W\, P_{eq} = P_{eq}\, , ~~{\rm in\ components}~~
 P_{eq}^{(l)}=\sum_k W^{(l)(k)} P_{eq}^{(k)}\ .
\end{equation}
{\r Theorem:} Under conditions (i), (ii) and (iii) the Boltzmann 
ensemble is the {\r only} equilibrium ensemble of the Markov process
and an {\r attractive fixed-point}.

\newpage

{\r Proof:} Let us first  define a distance between ensembles. Suppose 
we have two ensembles $\s E$ and $\s E'$, each of which is a collection 
of many configurations. Denote the probability for configuration $\s k$ 
in $\s E$ by $\s P^{(k)}$ and in $\s E'$ by $\s {P'}^{(k)}$. We define 
the {\r distance} between $\s E$ and $\s E'$ to be
\begin{equation} \label{distance} \s
 || E - E' ||\ =\ \sum_k | P^{(k)} - {P'}^{(k)} |\ , 
\end{equation}
where the sum goes over all configurations. Suppose that $\s E'$
resulted from the application of the transition matrix $\s W$ to the
ensemble $\s E$. We can compare the distance of $\s E'$ from the 
Boltzmann ensemble with the distance of $\s E$ from the Boltzmann 
ensemble: \s
\begin{eqnarray} 
  & &|| E' - E_B ||\ =\ \sum_l\, \left|\, \sum_k W^{(l)(k)}\, 
  (P^{(k)}-P_B^{(k)})\, \right|\ \qquad ({\rm using\ balance}) 
  \nonumber\\ &\le & 
  \sum_l \sum_k \left|\, W^{(l)(k)}\, (P^{(k)}-P_B^{(k)})\, \right|\ 
  \qquad ({\rm using\ the\ triangle\ inequality}) \nonumber\\ \s
  &=& \sum_k \left|\, P^{(k)}-P_B^{(k)} \,\right|\qquad 
  ({\rm using}\ W^{(l)(k)}\ge 0\ {\rm and\ normalization})\nonumber \\ 
  &=& ||E-E_B||\ . \label{contraction} \s
\end{eqnarray} \b
The last line is obtained by making use of the condition $\s 
\sum_l W^{(l)(k)}=1$ \b and of $\s W^{(l)(k)}\ge 0$. This shows that 
the Markov process can reduce the distance between an ensemble $\s E$
and the Boltzmann ensemble $\s E_B$ and will never increase this 
distance. Ergodicity is needed to proof that the Boltzmann ensemble 
is the only equilibrium ensemble and to rule out the equality in the 
$\s \le$ sign in equation (\ref{contraction}). We now study the 
approach to the fixed point and derive and explicit equation for 
the convergence.

The matrix $\s W$ by itself is normally not ergodic. Instead, due to 
the finiteness of the system, a number $\s n$ exists, so that the matrix 
$\s {\cal W}=W^{n}$ is ergodic. This means, all matrix elements of 
$\s {\cal W}$ are larger than zero. In particular, excluding states 
with infinite energy from our considerations,

\newpage

\begin{equation} \label{w_min} \s
  1 > w_{\min} = \min_{k,l} \left({\cal W}^{(k)(l)} \right) > 0
\end{equation}
holds. Let us assume that the state vectors $\s P$ and $\s P'$ are 
related by
\begin{equation} \s
  P'\ =\ {\cal W}\, P\,.
\end{equation}
As shown in the following, this implies the inequality
\begin{equation} \r
  || E' - E_B ||\ \le\ (1-w_{\min})\, || E - E_B ||\ .
\end{equation}
Let $\s\epsilon=||E-E_B||$. We can decompose the contributions to
$\s ||E-E_B||$ in $\s P^{(k)}-P_B^{(k)}\ge 0$ and in $\s
P^{(k)}-P_B^{(k)}<0$ terms,
\begin{equation} \label{decompose} \s
  \epsilon = ||E-E_B||\ 
             =\ \sum_{k\in K^+}\left(P^{(k)}-P_B^{(k)}\right) 
             +  \sum_{k\in K^-}\left(P_B^{(k)}-P^{(k)}\right)\,.
\end{equation}
Then the normalization $\s\sum_k P^{(k)}=\sum_k P_B^{(k)}=1$ implies
$$\s \epsilon - 2 \sum_{k\in K^-}\left(P_B^{(k)}-P^{(k)}\right) =
     \sum_{k\in K^+}\left(P^{(k)}-P_B^{(k)}\right) + 
     \sum_{k\in K^+}\left(P_B^{(k)}-P^{(k)}\right) = 0         $$
and similar for $\s \epsilon - 2 \sum_{k\in K^+}\left(P^{(k)} -
P_B^{(k)}\right).$ Therefore,
\begin{equation} \label{Psum_epsilon}\s
  \sum_{k\in K^+}\left(P^{(k)}-P_B^{(k)}\right) =
  \sum_{k\in K^-}\left(P_B^{(k)}-P^{(k)}\right) = \epsilon/2
\end{equation}
Using $\s{\cal W} P_B=P_B$, we have
(back in original order to apply $\s{\cal W}$)
{ \s \begin{eqnarray} \label{W_decomposition} 
&~& \epsilon = ||E'-E_B||\ = \\ \nonumber &~& \sum_l \left|\,
  \sum_{k\in K^+} {\cal W}^{(l)(k)}\,\left( P^{(k)}-P_B^{(k)}\right) -
  \sum_{k\in K^-} {\cal W}^{(l)(k)}\,\left( P_B^{(k)}-P^{(k)}\right)\,
  \right| \\ \nonumber \s
&=& \sum_{l\in L^+} \left( 
    \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)}\right) 
  - \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right) 
    \right) \\ \nonumber \s
&+& \sum_{l\in L^-} \left( 
    \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right) 
  - \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)} \right)
    \right)
\end{eqnarray} }
Here $\s L^+$ is defined as the set of configurations $\s l$ which 
fulfills the inequality
$$\s \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)}\right) 
 \ge \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right)\ ,$$
while $\s L^-$ is the set of configurations $\s l$ with
$$\s \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right)
 > \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)}\right)\ .$$
Normalization (\ref{normalization}) and that the smallest matrix 
element is larger than zero (\ref{w_min}) imply
$$\s 1\ \ge\ \sum_{l\in L^+} {\cal W}^{(l)(k)}\ \ge\ w_{\min} $$
and the same equation with $\s L^+$ replaced by $\s L^-$. Inserting this
into equation (\ref{W_decomposition}), picking the extremes $\s 1$ 
and $\s w_{\min}$ appropriately, yields 
\begin{equation} \label{W_inequality} \s
  ||E'-E_B|| \le (1-w_{\min})\,\epsilon/2\ 
  +\ (1-w_{\min})\,\epsilon/2\ =\ (1-w_{\min})\,\epsilon\ .
\end{equation}
Under repeated application of the matrix we obtain a state
\begin{equation} \label{MC_chain} \s
 P^t\ =\ {\cal W}^t\,P,~~{\rm with}~~ t=1,2,\dots 
\end{equation}
and find for the {\r approach to the equilibrium ensemble}
\begin{equation} \label{MC_convergence} \s
  ||E^t - E_B||\ \le\  \exp ( - \lambda\,t)\,||E-E_B||\,,~~
  \lambda=-\ln(1-w_{\min})>0\,.
\end{equation}
{\r Hence, the equilibrium ensemble is an attractive fixed-point.}
\vfill \newpage

There are many ways to construct a Markov process satisfying (i), (ii)
and (iii).  A stronger condition than balance~(\ref{balance}) is
%
\begin{description}
  \item{(iii')} {\r Detailed balance}:
\begin{equation} \label{detailed_balance} \s
  W^{(l)(k)}\, e^{-\beta E^{(k)}}\ =\ W^{(k)(l)} e^{-\beta E^{(l)}}\ .
\end{equation}
Using the normalization $\s \sum_k W^{(k)(l)} = 1$ detailed
balance implies balance (iii).
\end{description}
%
At this point we have replaced the canonical ensemble average by a time 
average over an artificial dynamics. Calculating then averages over large 
times, like one does in real experiments, is equivalent to calculating 
ensemble averages. One distinguishes {\r dynamical universality classes}.  
The Metropolis and heat bath algorithms discussed in the following fall 
into the class of {\r model~A or Glauber dynamics}, which imitates the 
thermal fluctuations of nature to some extent. Cluster algorithms 
discussed constitute another universality class. Some recent attention 
has focused on dynamical universality classes of non-equilibrium systems.

\foilhead{\r The Metropolis Algorithm} \vskip -30pt

Detailed balance still does not uniquely fix the transition
probabilities $\s W^{(l)(k)}$. The Metropolis algorithm can 
be used whenever one knows how to calculate the energy of a 
configuration. Given a configuration $\s k$, the Metropolis 
algorithm proposes a configuration $\s l$ with probability
\begin{equation} \label{Metropolis_f} \s
 f(l,k) ~~{\rm normalized\ to}~~~ \sum_l f(l,k)=1\ .
\end{equation}
For $\s f(l,k)$ we derive a symmetry condition which ensures detailed 
balance.  The new configuration $\s l$ is accepted with probability
\begin{equation} \label{Metropolis} \s
  w^{(l)(k)} = \min \left[ 1,\, {P_B^{(l)}\over P_B^{(k)}} \right] 
  = \left\{ \matrix{ 1& {\rm for}& E^{(l)}<E^{(k)}\cr
  e^{-\beta (E^{(l)}-E^{(k)})} &{\rm for}& E^{(l)}>E^{(k)}.\cr}\right.
\end{equation}
If the new configuration is rejected, the {\r old configuration has to 
be counted again.}

The Metropolis procedure gives rise to the transition probabilities \s
\begin{eqnarray} \label{Metropolis_W_ne}
  W^{(l)(k)} &=& f(l,k)\, w^{(l)(k)}~~{\rm for}~~l \ne k\\ {\rm and}~~
  W^{(k)(k)}& =& f(k,k) + \sum_{l\ne k} f(l,k)\, (1-w^{(l)(k)})\ .
\end{eqnarray} \b
Therefore, the ratio $\s \left( W^{(l)(k)}/W^{(k)(l)} \right)$ 
satisfies detailed balance~(\ref{detailed_balance}) if
\begin{equation} \label{symmetry} \s
  f(l,k)\ =\ f(k,l) ~~{\rm holds}\,.
\end{equation}
Otherwise the probability density $\s f(l,k)$ is unconstrained. So there 
is an amazing flexibility in the choice of the transition probabilities 
$\s W^{(l)(k)}$.  One can even use acceptance probabilities distinct 
from those of equation~(\ref{Metropolis}) and the proposal probabilities 
are then not necessarily symmetric anymore {\r (Hastings}).} Also, 
the algorithm generalizes immediately to arbitrary weights.

\noindent
The {\r acceptance rate} is defined as the ratio of accepted {\r 
changes} over proposed moves (moves proposing the at hand 
configuration are not counted as accepted).

\newpage

{\r\bf Potts Model Heatbath Algorithm (Gibbs Sampler)}

The {heatbath algorithm} chooses a state $\s q_i$ directly with the 
local Boltzmann distribution defined by its nearest neighbors. The 
state $\s q_i$ can take on one of the values $\s 1,\dots,q$ and, with 
all other states set, determines a value of the energy function. We 
denote this energy by $\s E(q_i)$ and the Boltzmann probabilities are
\begin{equation} \label{Potts_B_weights} \s
  P_B(q_i)\ =\ {\rm const}\ e^{-\beta\, E(q_i)}
\end{equation}
where the constant is determined by the normalization condition
\begin{equation} \label{Potts_B_normalization} \s
  \sum_{q_i=1}^q P_B(q_i)\ =\ 1\ .
\end{equation}        
In equation~(\ref{Potts_B_weights}) we can define $\s E(q_i)$ to be just 
the contribution of the interaction of $\s q_i$ with its nearest neighbors 
to the total energy and absorb the other contribution into the overall 
constant. The $\s E(q_i)$ values depend only on how the nearest neighbors 
of the spin $\s q_i$ partition into the values $\s 1,\dots,q$. For low values 
of $\s q$ and the dimension $\s d$ the most efficient implementation of the 
heatbath algorithm is to tabulate all possibilities. However, here we 
prefer to give a generic code which works for arbitrary values of $\s q$ 
and $\s d$.

For this we calculate the cumulative distribution function of the 
heat bath probabilities
\begin{equation} \label{prob_hb} \s
  P_{HB}(q_i) = \sum_{q'_i=1}^{\tt q_i} P_B(q'_i)\ .
\end{equation}
The normalization condition~(\ref{Potts_B_normalization}) implies
$\s P_{HB}(q)=1$. Comparison of these cumulative probabilities
with a uniform random number $\s x^r$ yields the heat bath update
$\s q_i\to q'_i$. Note that in the heat bath procedure the original value
$\s q_i^{\rm in}$ does not influence the selection of $\s q_i^{\rm new}$. 

\newpage

{\r\large\bf Start and Equilibration}

Initially we have 
to start with a microstate which may be far off the Boltzmann 
distribution. Far off means, that the Boltzmann probability (at
temperature $\s T$) for a typical state of the initially generated
distribution can be very, very small. Suppression factors like
$\s 10^{-10000}$ are well possible. Although the weight of states 
decreases with $\s 1/n$ where $\s n$ is the steps of the Markov process,
one should exclude the initial states from the equilibrium statistics. 
In practice this means we should allow for a certain number of sweeps 
{\tt nequi} to equilibrate the system.

Many ways to generate start configurations exist. Two natural and
easy to implement choices are:

\begin{enumerate}
\item Generate a random configuration corresponding to $\s \beta =0$.
This defines a {\r random or disordered start} of a MC simulation.

\item Generate a configuration for which all Potts spins take on the 
same $\s q$-value. This is called an {\r ordered start} of a MC simulation.
\end{enumerate}

\newpage
{\bf Examples of initial time series:}

 \begin{figure} \begin{center}
 \includegraphics[width=0.7\textheight]{figures/2dI_ts.pdf}
 \includegraphics[width=0.7\textheight]{figures/2d10q_ts.pdf}
 \caption{\small Left: Two Metropolis time series of 200 sweeps each 
for a $\s 2d$ Ising model on a $\s 80\times 80$ lattice at $\s \beta=0.4$ 
are shown. Ordered and disordered starts are used. The exact mean value 
$\s \widehat{e}_{0s}=-1.10608$ is also indicated (assignment
{\tt a0302\_01}). Right: $\s q=10$ Potts model time series of 200 
sweeps on an $\s 80\times 80$ lattice at $\s \beta=0.62$. Measurements 
of the action variable after every sweep are plotted 
for ordered and disordered starts (assignment {\tt a03003\_05}). }
 \end{center} \end{figure}

\foilhead{\r Computer Program } \vskip -30pt

\tt\s

program potts\_ts.f

~~~~~subroutine potts\_init.f

~~~~~~~~~~subroutine lat\_init.f

~~~~~~~~~~subroutine potts\_act\_tab.f 

~~~~~~~~~~subroutine potts\_order.f 

~~~~~~~~~~subroutine potts\_ran.f 

~~~~~~~~~~subroutine potts\_act.f 

~~~~~~~~~~subroutine potts\_wght.f 

~~~~~subroutine potts\_met.f \rm

\foilhead{\r Energy Checks} \vskip -30pt \b

For the {\r $\s 2d$ Ising model} we can test against the exact finite 
lattice results of {\r Ferdinand and Fisher.} We simulate an $\s 20^2$ 
lattice at $\s \beta=0.4$. 

Now we use a statistics of \s 10$\,$000 sweeps for reaching equilibrium 
\b and assume that this is an overkill. A more careful analysis is
the subject of the next lecture.
 
The statistics for \s measurement \b is chosen 32 times larger: \s 
64 bins of 5$\,$000 sweeps \b each. The number 64 is taken, because 
according to the student distribution the approximation to the Gaussian 
approximation is then already excellent, while the binsize of 5$\,$000 
($\gg 200$) is argued to be large enough to neglect correlations 
between the bins.  With this statistics we find (assignment 
{\tt a0303\_06})
\begin{equation} \label{2dI_es} \s
\overline{e}_{0s} = -1.1172\ (14)~~ {\rm (Metropolis)} 
~~~{\rm versus}~~~ \widehat{e}_s = -1.117834~~ {\rm (exact)}\ .
\end{equation}
Performing the Gaussian difference test gives a perfectly admissible 
value ${\r Q=0.66 .}$ 

\foilhead{\r Specific Heat } \vskip -30pt

With $\s \widehat{E}=\langle E\rangle$ the
{\bf specific heat} is defined by
$$\s C\ =\ {d\,\widehat{E}\over d\,T}\ =\ \beta^2\, \left( \,
     \langle E^2 \rangle - \langle E \rangle^2 \right)\ .$$
The last equal sign follows by working out the temperature derivative 
in the definition of the mean energy
and is known as {\bf fluctuation-dissipation theorem}. It is often 
used to estimate the specific heat from equilibrium simulations without 
relying on numerical derivatives. An equivalent formulation is
$$\s C\ =\ \beta^2\, \left\langle \left( E - \widehat{E}\right)^2
                      \right\rangle $$
as is easily shown by expanding $\s \left(E-\widehat{E}\right)^2$ and 
working out the expectation values. Defining
$$\s {\tt act2l}=\overline{\left(\frac{\tt iact}{\tt mlink}\right)^2}$$
we have a notation close to the computer code:
$$\s C\ =\ \frac{\beta^2\,N\,d^2}{n}\,\sum_{i=1}^n
     \left({\tt act2l}_i-{\tt actlm}^2\right)\,,$$
or
$$\s C\ =\ \frac{\beta^2\,N\,d^2}{n}\,
     \sum_{i=1}^n ({\tt actl}_i-{\tt actlm})^2\,,$$
where the sums are over all measurements in the times series. 
When energy {\r histograms} are available, this is calculated 
by equations in {\tt\s potts\_mu2.f} of {\tt ForLib}.

In the limit of an infinite statistics specific heat estimates
from all these equation agree. But with a finite statistics a number 
of problems emerge. In the binning approach, where {\s\tt nrpt} of 
the production run defines the number of blocks, one may want to 
use for $\s \widehat{E}$ estimators $\s \overline{E}_i$ which are 
constructed from the histograms of the blocks. With $\s N_b$ the 
number of data in each bin and $\s H_i$ the energy histogram of 
block $\s i$:
$$\s \overline{E}_i = \frac{1}{N_b} \sum_{j \in {\rm block}(i)} E_j
     \ = \frac{\sum_E E\,H_i(E)}{\sum_E H_i(E)}\ .$$
Estimates of $\s \overline{C}_i$ from either equation agree then and 
Gaussian error bars are expected due to the binning. However, for a 
not so good statistics a bias towards too small $\s \overline{C}_i$ 
values occurs, because $\s E_j$ and $\s \overline{E}_i$ come from 
the same block in the estimate
$$\s  \overline{C}_i = \frac{\beta^2}{N_b} \sum_{j\in{\rm block}(i)} 
  \left( E_j^{\,2} - \overline{E}_i^{\,2} \right) =  
  \frac{\beta^2\, \sum_E \left(E^{\,2}-\overline{E}_i^{\,2}\right)
  \,H_i(E)}{\sum_E H_i(E)} $$
$$\s   = \frac{\beta^2}{N_b} \sum_{j\in{\rm block}(i)} 
  \left(E_j-\overline{E}_i\right)^2 = \frac{\beta^2\,
  \sum_E (E-\overline{E}_i)^2\,H_i(E)}{\sum_E H_i(E)}\ . $$
The $\s \overline{E}_i$ estimators are certainly inferior to the estimate 
$\s \overline{E}$, which relies on the entire statistics. Therefore, one
may consider to use $\s \overline{E}$ instead of $\s \overline{E}_i$. 
However, then one does not know anymore how to calculate the error bar 
of $\s \overline{C}$ as the $\s \overline{C}_i$ estimates would rely 
on overlapping data. The situation gets even worse when we include 
reweighting, as this non-linear procedure implies that the two 
estimators will in general differ. These difficulties are overcome 
by the {\r jackknife method.}

When histograms are used the fast way to create jackknife bins is 
to sum first the entire statistics:
$$\s  H(E) = \sum_{i=1}^{N_b} H_i(E)\ .$$
Subsequently jackknife histograms (superscript $\s J$) are defined by
$$\s  H_i^J(E) = H(E) - H_i(E)$$
and jackknife estimates $\s \overline{C}_i^J$ are obtained by using 
$\s H_i^J(E)$ instead of $\s H_i(E)$.

\vfil
\end{document}

