
% BB Aug 30, 2006.
% BB Aug 28, 2005.
% BB May 24, 2005.
% BB Dec 16, 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} Leipzig University, December 16, 2004}
\date{\color{blue} FSU,  August 31, 2006}
\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}}

$$ u(x) = \cases{ 1~~{\rm for}~~0\le x < 1;\cr
                  0 ~~{\rm elsewhere}. } $$ \b
The corresponding distribution function is \s
$$ U(x) = \int_{-\infty}^x u(x')\,dx' = \cases{
   0~~{\rm for}~~x<0;\cr
   x~~{\rm for}~~0\le x\le 1;\cr 1~~{\rm for}~~x>1.} $$ \b
It allows for the construction of general probability distributions.  
Let \s
$$ y = F(x) = \int_{-\infty}^x f(x')\,dx'\ . $$ \b
For \s $y^r$ \b being a uniformly distributed random variable in \s $[0,1)$ 
$$ x^r = F^{-1} (y^r)~~{\rm\b is\ then\ distributed\ according\ to\
   the\ probability\ density}~~f(x)\ . $$

\foilhead{\b Example } \b \vskip -20.0pt

\noindent \centerline{Mapping of the uniform to the Cauchy distribution.}

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/gdf.pdf}
% \caption{\small Mapping of the uniform to the Cauchy distribution.}
\end{center} \end{figure}

\foilhead{\r Confidence Intervals and Sorting}\b \vskip -20pt

\noindent One defines {\bf ${\bf q}$-tiles} (also {\bf quantiles} or 
{\bf fractiles}) $\s x_q$ of a distribution function by 
$$\s F(x_q)\ =\ q\ . $$
An example is the {\bf median} $\s x_{1\over 2}.$ The probability 
content of the {\bf confidence interval}
$$\s [x_q,x_{1-q}] ~~{\rm is}~~ p=1-2q\ .$$
Example: {\r Gaussian or normal distribution} of variance $\sigma^2:$
$$\s [-n\sigma,+n\sigma]\ \Rightarrow\ p=0.6827\ {\rm for}\ n=1,~~
     p=0.9545\ {\rm for}\ n=2\ .$$
The {\bf peaked distribution function} 
$$\s F_q(x) = \cases{ F(x)\ {\rm for}\ F(x) \le {1\over 2},\cr
                 1 - F(x)\ {\rm for}\ F(x) > {1\over 2}.} $$
provides a graphical visualization of probability intervals of
such a distribution:
%
\begin{figure} \begin{center}
\includegraphics[height=0.5\textheight]{figures/gqdf.pdf}
% \caption{\small Gaussian peaked distribution function and estimates 
% of $x_q$ for the 70\% (approximately $1\,\sigma$) and 95\% 
% (approximately $2\,\sigma$) confidence intervals.}
\end{center} \end{figure}

\noindent {\r Sorting} allows for an empirical estimate. Assume we 
generate $n$ random number $\s x_1,...,x_n$. We may re-arrange the 
$\s x_i$ in increasing order:
$$\s x_{\pi_1} \le x_{\pi_2} \le \dots \le x_{\pi_n} $$
where $\s \pi_1,\dots , \pi_n$ is a permutation of $\s 1,\dots ,n$.

\noindent An estimator for the distribution function $\s F(x)$ is 
then the {\bf empirical distribution function}
$$\s \overline{F} (x) = {i\over n}~~~{\rm for}~~~x_{\pi_i} \le x <
     x_{\pi_{i+1}},\ i=0, 1,\dots , n-1, n $$
with the definitions $\s x_{\pi_0}=-\infty$ and 
$\s x_{\pi_{n+1}}=+\infty$. To calculate $\overline{F}(x)$ one needs 
an efficient way to {\bf sort} $n$ data values in ascending (or 
descending) order. In the STMC package this is provided by a 
{\bf heapsort} routine, which arrives at the results in
$\s O(n\,\log_2 n)$ steps.

\noindent Example: {\r Gaussian distribution in assignment 
          {\tt a0106\_02}} (200 and 20,000 data). \vskip 20pt

\newpage

{\large\r Central Limit Theorem:}  Convergence of the Sample Mean
$$\s  {\rm Gaussian}~~~\sigma^2 (\overline{x}^r) = 
      {\sigma^2 (x^r)\over N}~~~{\rm for}~~~
      \overline{x} = \frac{1}{N} \sum_{i=1}^N x_i\ .$$

\foilhead{\r Binning } \b \vskip -20pt

\noindent We group {\tt NDAT} data into {\tt NBINS} bins, where each 
binned data point is the arithmetic average of
$$\s {\tt NBIN=[NDAT/NBINS]}~~({\rm Fortran\ integer\ division.}) $$
original data points.  Should not be confused with histogramming! 
Preferably {\tt NDAT} is a multiple of {\tt NBINS}. The purpose of 
the binning procedure is twofold:

\begin{enumerate}
\item When the the central limit theorem applies, the binned data
will become practically Gaussian when  {\s\tt NBIN} is large enough. 
This allows to apply \r Gaussian error analysis \b methods even when 
the original are not Gaussian.

\item When data are generated by a Markov process subsequent events 
are correlated. For binned data these \r autocorrelations are reduced 
\b and can be neglected, once {\s\tt NBIN} is sufficiently large.
\end{enumerate}

\newpage
{\bf Example:}

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/bin_his.pdf}
\caption{\small Comparison of a histogram of 500 binned data with 
the normal distribution $\sqrt{(120/\pi)}\, \exp [-120\,(x-1/2)^2]$.
Each binned data point is the average of 20 uniformly distributed 
random numbers. \r Assignment {\tt a0108\_02}.}
\end{center} \end{figure}

\foilhead{\r Gaussian Difference Test} \b \vskip -29.0pt

\noindent One is often faced with comparing two different empirical 
estimates of some mean. How large must 
$\s D={\overline x}-{\overline y}$ be in order to indicate a 
real difference? The quotient 
$$\s d^r\ =\ {D^r/\sigma_D}\,,\qquad \sigma_D^2=
     \sigma_{\overline{x}}^2+\sigma_{\overline{y}}^2   $$
is normally distributed with expectation zero and variance one.
The {\bf likelihood that the observed difference 
$\s |{\overline x}-{\overline y}|$ is due to chance}
is defined to be
$$\s Q\ =\ 1-P(|d^r|\le d)\ =\ 2\,G_0(-d)\ =\ 
1-{\rm erf} \left( {d/\sqrt{2}}\right)\ .$$
If the assumption is correct, then $\s Q$ is a 
uniformly distributed random variable in the range $\s [0,1)$.
{\bf Examples:} {\r (interactive in {\tt ForProc/Gau\_dif/})}
%
\begin{table}[ht]
% \caption{Gaussian difference tests (assignment {\tt a0201\_06}).} 
\medskip
\centering
\begin{tabular}{||c|c|c|c|c|c||}                            \hline
 ${\overline x}_1 \pm \sigma_{{\overline x}_1}$ & $1.0\pm 0.1$
 & $1.0\pm 0.1$ &$1.0\pm 0.1$&$1.0\pm 0.05$ & $1.000\pm 0.025$
\\ \hline
 ${\overline x_2} \pm \sigma_{{\overline x}_2}$ & $1.2\pm 0.2$
 &$1.2\pm 0.1$ &$1.2\pm 0.0$&$1.2\pm 0.00$ & $1.200\pm 0.025$
\\ \hline
 $Q$ & 0.37 & 0.16 &~~~0.046 & 0.000063 & $0.15\times 10^{-7}$
\\ \hline
\end{tabular}
\label{tab_gdt}
\end{table}

\newpage
\noindent
{\r\large\bf Gaussian Error Analysis for Small Samples:} \b

\noindent
Gosset's Student Distribution and the Student Difference Test.

\noindent
{\r\large\bf The Error of the Error Bar:} \b

\noindent
$\s\chi^2$ Distribution and the Variance Ratio Test (F-Test).

\noindent {\r\large\bf The Jackknife Approach:} 

\noindent
Jackknife estimators correct for the bias and the error of the bias. 
If there is not bias, their results are identical with the conventional
error analysis. As the extra effort of using Jacknife routines is
minimal, they should be the {\r standard of error analysis.}

\noindent
Bias problems occur when one estimates a non-linear function of a 
mean \s ${\widehat x}$: 
$$ {\widehat f} = f({\widehat x})\ . $$
\b Jacknife estimators of the function are then defined by 
$$\s {\overline f^J}\ =\ {1\over N} \sum_{i=1}^N f^J_i
   ~~{\rm with}~~ f_i^J\ =\ f(x_i^J) ~~{\rm and}~~ 
   x_i^J\ =\ {1\over N-1} \sum_{k\ne i} x_k\ .  $$
The estimator for the variance is
$$\s s^2_J ({\overline f^J})\ =\ {N-1 \over N} 
     \sum_{i=1}^N (f^J_i - {\overline f^J})^2\ . $$
Straightforward algebra shows that in the unbiased case the jackknife 
variance reduces to the normal variance. Of order 
$\s N$ (not $\s N^2$) operations are needed to construct the jackknife 
averages $\s x_i^J,\ i=1,\dots ,N$ from the original data.

\noindent
{\r\large\bf Determination of Parameters (Fitting):} \b

\noindent Levenberg-Marquardt approach besides simple linear
regression.


\foilhead{\r Statistical Physics and Markov Chain Monte Carlo 
             Simulations} \vskip -30pt

\noindent \b MC simulations of systems described by the Gibbs canonical 
ensemble aim at calculating estimators of physical observables at a 
temperature $\s T$. In the following we consider the calculation of 
the {\r expectation value} of an {\r observable} $\s \mathcal{O}$. All 
systems on a computer are discrete, because a finite word length has 
to be used. Hence, 
{\s \begin{eqnarray} \nonumber
  \widehat{\mathcal{O}} = \widehat{\mathcal{O}} (\beta) = 
  \langle \mathcal{O} \rangle &=& Z^{-1} \sum_{k=1}^K
  \mathcal{O}^{(k)}\,e^{-\beta\,E^{(k)} }\\ \nonumber
  {\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$. 
% \noindent {\r Continuous systems} are also treated in the book.

\newpage
In the following I use {\r Potts models} on $d$-dimensional cubic 
lattices with periodic boundary conditions. Without being overly 
complicated, these models allow to illustrate essential features 
we are interested in. We define the energy of the system by
$$\s E_0^{(k)} = - 2\, \sum_{<ij>} \delta (q_i^{(k)},q_j^{(k)}) + 
  {2\,d\,N\over q}~~~{\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.\, .  $$ 
The sum $\s <ij>$ is over the nearest neighbor lattice sites and 
the Potts {\r spins} or {\r states} of $\s q_i^{(k)}$ take the values 
$\s 1,\dots ,q$. For the {\r energy per spin} the notation is
$\s e_s = E/N\,$ and our normalization is chosen so that $\s e_s$ 
agrees for $\s q=2$ with the conventional Ising model definition.
For the $\s 2d$ Potts models a number of exact results are known,
e.g.,
$$ \s \beta_c\ =\ {1\over T_c} = {1\over 2}\, \ln (1 + \sqrt{q} )
   = \beta_c^{\rm Potts}\,,\qquad
  e_{0s}^c = E_0^c/N = {4\over q} - 2 - 2/\sqrt{q}\ . $$
% = {4\over q} + e_{0s}[{\rm Baxter}]
The phase transition is second order for $\s q\le 4$ and first order 
for $\s q\ge 5$.

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

\noindent A Markov chain allows to generate configurations $\s k$
with probability
$$\s  P_B^{(k)} = c_B\, w^{(k)}_B = c_B\, e^{-\beta E^{(k)}}\,,
   \qquad c_B~~{\rm constant}\,. $$
The {\r state vector} $\s (P_B^{(k)})$, for which the configurations are 
the vector indices, is called {\r Boltzmann state}. A Markov chain is a 
simple dynamic process, which generates configuration $\s k_{n+1}$
stochastically from configuration $\s k_n$. 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]$. Then, the transition matrix
$$\s W = \left( W^{(l)(k)}\right) $$
defines the Markov process. Note, that this matrix is a very big 
and never stored in the computer. The matrix achieves our goal and 
generates configurations with the desired probabilities, when it
satisfies the following properties:
\newpage

\begin{description}

\item{(i)} {\r\bf Ergodicity}:
$$\s  e^{-\beta E^{(k)}} > 0 ~~{\rm and}~~ e^{-\beta E^{(l)}} > 0
  ~~~{\rm imply:} $$
an integer number $\s n>0$ exists so that $\s (W^n)^{(l)(k)}>0$ holds.

\item{(ii)} {\r\bf Normalization}:
$$\s \sum_l W^{(l)(k)} = 1\ . $$

\item{(iii)} {\r\bf Balance}: 
$$\s \sum_k W^{(l)(k)}\,e^{-\beta E^{(k)}}\ =\ e^{-\beta E^{(l)}}\ .$$
Balance means: The Boltzmann state is an eigenvector with eigenvalue 
$\s 1$ of the matrix $\s W = (W^{(l)(k)})$.

\end{description}
% \newpage

\noindent With that we have replaced the canonical ensemble average by 
a time average over an artificial dynamics and one
distinguishes {\it dynamical universality classes}.  

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

\noindent Detailed balance still does not uniquely fix the transition
probabilities $\s W^{(l)(k)}$. The Metropolis algorithm is a popular 
choice 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
$$\s f(l,k)=f(k,l)~~{\rm normalized\ to}~~\sum_l f(l,k)=1\ .$$
The new configuration $l$ is accepted with probability
$$\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.$$
If the new configuration is rejected, the old configuration has to 
be counted again.  % A beginners mistake is to count configurations 
%                    only after accepted changes.      \newpage

\foilhead{\r The Heatbath algorithm} \vskip -30pt \b

The heat bath algorithm chooses a spin $q_i$ directly with the 
local Boltzmann distribution defined by its nearest neighbors 
$$\s P_B(q_i)\ =\ {\rm const}\ e^{-\beta\, E(q_i)}\ . $$
As many Metropolis hits (on the same spin) are needed to reach
this distribution, the heatbath is more efficient than the 
Metropolis algorithm. However, it models more complicated than
Potts models the calculation of the local heatbath probabilities 
is often too involved to make it a viable alternative.
\bigskip
\bigskip

\centerline{\r\large\bf Start and equilibration} \b
\noindent Initially we have to start with a microstate which may be 
far off the Boltzmann distribution. 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. 
Many ways to generate start configurations exist, e.g.,

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

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

{\bf Examples:} (assignments {\tt a0303\_01} and {\tt a0303\_05})

\begin{figure} \begin{center}
\includegraphics[width=0.75\textheight]{figures/2dI_ts.pdf}
\includegraphics[width=0.75\textheight]{figures/2d10q_ts.pdf}
\caption{\small Initial time series of 200 sweeps each on a 
$80\times 80$ lattice. Left: Metropolis for a $2d$ Ising model 
at $\beta=0.4$.
% The exact mean value $\widehat{e}_{0s}=-1.10608$ is also indicated 
% (assignment {\tt a0302\_01}). 
Right: $q=10$ Potts model at $\beta=0.62$. }
% (assignment {\tt a03003\_05}). 
\end{center} \end{figure}

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

For the {\r $2d$ Ising model} we can test against the exact finite 
lattice results of {\r Ferdinand and Fisher.} We simulate a $\s 20^2$ 
lattice at $\s \beta=0.4$ using  {\s 10$\,$000} sweeps for reaching 
equilibrium {\s 64 bins of 5$\,$000 sweeps} for measurement (a careful 
justification is given later). We find {\r (assignment {\tt a0303\_06})}
$$\s \overline{e}_{0s} = -1.1172\ (14)~~ {\rm (Metropolis)} 
~~~{\rm versus}~~~\widehat{e}_s = -1.117834~~ {\rm (exact)}\ .$$
The {\r Gaussian difference test} gives a perfectly admissible value 
$\r Q=0.65\ .$ 

For the {\r $2d$ 10-state Potts model} at $\s \beta =0.62$ we test our 
Metropolis versus our heat bath code on a $20\times 20$ lattice to
find (assignment {\tt a0303\_08})
$$\s  {\tt actm} = 0.321772\ (75)~~{\rm (Metropolis)}~~{\rm versus}
    ~~~{\tt actm} = 0.321661\ (70) ~~{\rm (heat\ bath)} $$
and $\s Q=0.28$ for the Gaussian difference test. Another perfectly
admissible value.

\foilhead{\r $3d$ 3-state Potts model} \vskip -30pt \b
\noindent
To illustrate features of a first order phase transition we 
simulate the $3d$ 3-state Potts model on a $\s 24^3$ lattice at
a {\r pseudo-transition temperature} and plot its energy
histogram. A characteristic {\r double peak} structure is 
found (assignment {\tt a0303\_10}):
%
\begin{figure} \begin{center}
\includegraphics[height=0.55\textheight]{figures/3d3q_h.pdf}
\end{center} \end{figure}
\newpage
%
\newpage

{\r Self-Averaging Illustration for the $2d$ Heisenberg model}

We compare the peaked distribution function of a mean energy density
per link for different lattice sizes. The property of {\bf self-averaging} 
is observed: The larger the lattice, the smaller the confidence range.
The other way round, the peaked distribution function is very well 
suited to exhibit observables for which self-averaging does not work, 
as for instance encountered in spin glass simulations.
%
\begin{figure} \begin{center}
\includegraphics[height=0.45\textheight]{figures/qdf_o3.pdf}
\caption{$O(3)$ $\sigma$-model at $\beta=1.1$ (assignments 
{\tt a0304\_06} and {\tt a0304\_08}).}
\end{center} \end{figure}

\foilhead{\r Statistical Errors of Markov Chain MC Data}\vskip -30pt

\b \noindent A typical MC simulation falls into two parts: 

\begin{enumerate}

\item Equilibration without measurements.

\item Production with measurements.

\end{enumerate}

\noindent Rule of thumb (for long calculations): {\bf Do not spend more than 
50\% of your CPU time on measurements!} 
\hfil\break

\centerline{\r\large\bf Autocorrelations}
% \foilhead{\r Autocorrelations} \vskip -30pt

\noindent 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.

\newpage 
\noindent Given is a {\r time series} of measurements $f_i$, 
$i=1,\dots,N$.  With the notation $\s t=|i-j|$ the autocorrelation 
function of the mean $\s \widehat{f}$ is defined by
$$\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} $$
Some algebra shows that the variance of the estimator 
$\s \overline{f}$ for the mean and the 
autocorrelation functions are related by
$$\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) = {\widehat{C}(t)\over \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} as
$$\s \tau_{\rm int}\ =\ \lim_{N\to\infty} \left[ 1+2\sum_{t=1}^{N-1} 
  \left(1-{t\over N}\right)\, \widehat{c}(t) \right]\ .$$
% \newpage

\centerline{\r\large\bf Self-consistent versus reasonable error analysis}
\bigskip % \vskip -30pt

% \noindent 
The calculation of the integrated autocorrelation provides a {\r 
self-consistent} error analysis. But in practice this is often 
out of reach. 

According to the Student distribution about twenty independent data 
are sufficient to estimate mean values reliably, while the error
of the error bar the $\chi^2$ distribution implies that about one
thousand are needed for an estimate of the integrated autocorrelation
time with {\s 10\% accuracy on the two $\sigma$ confidence level.}

In practice, one may rely on the binning method with a fixed number 
of {\s $\ge 16$ bins.} How do we know then that the statistics has 
become large enough? There can be indirect arguments like finite size 
scale extrapolations, which suggest that the integrated autocorrelation 
time is (much) smaller than the achieved bin length. This is no longer 
self-consistent, but a {\r reasonable error analysis}.

\foilhead{\r Comparison of Markov chain MC algorithms}
\vskip -30pt
\centerline{\bf The $d=2$ Ising model at the critical temperature}

\begin{figure} \begin{center}
\includegraphics[height=0.54\textheight]{figures/iat2_large.pdf}
\caption{One-hit Metropolis algorithm with sequential updating: 
{\r critical slowing down}, $\tau_{\rm int}\sim L^z$ where
$z\approx 2.17$ is the {\r dynamical critical exponent}
(assignment {\tt a0402\_02}~D). }
\end{center} \end{figure}
%

\newpage
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 
\r (assignment {\tt a0503\_05}). }
\end{center} \end{figure}

\foilhead{\r Simulations of the Multicanonical Ensemble} \b
\vskip -30pt

One of the questions which ought to be addressed before performing
a large scale computer simulation is {\r ``What are suitable weight
factors for the problem at hand?''} So far we used the Boltzmann 
weights as this appears natural for simulating the
Gibbs ensemble. However, a broader view of the issue
is appropriate.  

Conventional simulations can by re-weighting techniques only be 
extrapolated to a vicinity of the simulation temperature. For 
multicanonical simulations this is different. A single simulation 
allows to obtain equilibrium properties of the Gibbs ensemble
over a {\r range of temperatures.} Of particular interest are two 
situations for which canonical simulations do not provide the 
appropriate implementation of importance sampling: 
\begin{enumerate}
\item The physically important configurations are rare in the 
      canonical ensemble.
\item A rugged free energy landscape makes the physically important 
      configurations difficult to reach. 
\end{enumerate}

\newpage
Multicanonical simulations sample, in an appropriate energy range, with 
an {\r approximation} to the weights 
$$\s \widehat{w}_{mu}(k)=w_{mu}(E^{(k)}) = 
   e^{-b(E^{(k)})\, E^{(k)} + a (E^{(k)})} = {1\over n(E^{(k)})} $$
where $\s n(E)$ is the {spectral density}. The function $\s b(E)$ 
defines the inverse {\r microcanonical temperature} and $\s a(E)$ the 
{\r dimensionless, microcanonical free energy}. The function $\s b(E)$ 
has a relatively smooth dependence on its arguments, which makes it a 
useful quantity when dealing with the weight factors.  
The multicanonical method requires {\r two steps:}

\begin{enumerate}

\item Obtain a {\r working estimate} the weights. Working estimate means
      that the approximation has to be good enough so that the simulation
      covers the desired eneryg or temperature range. 

\item Perform a Markov chain MC simulation with the {\r fixed} weights.
      The thus generated configurations constitute the {\r multicanonical 
      ensemble}.

\end{enumerate}

\vfil
\newpage

\foilhead{\r Re-Weighting to the Canonical Ensemble}\vskip -30pt \b

Given the {multicanonical time series}, where $i=1,\dots,n$ labels 
the generated configurations. The formula 
$$\s \overline{\mathcal{O}} = { \sum_{i=1}^n \mathcal{O}^{(i)}\, \exp 
 \left[ -\beta\,E^{(i)} + b(E^{(i)})\,E^{(i)}-a(E^{(i)}) \right] \over 
 \sum_{i=1}^n \exp \left[ -\beta\,E^{(i)} + 
 b(E^{(i)})\,E^{(i)}-a(E^{(i)}) \right] }\ .$$
replaces the multicanonical weighting of the simulation by the Boltzmann 
factor. The denominator differs from the partition function $Z$ by a 
constant factor which drops out (for discrete systems this simplifies
for functions of the energy using histograms). 
The computer implementation of these equations requires care and a {\r 
Jackknife analysis} with {\r logarithmic coding} relying on the formula
$$\s \ln C = \max \left( \ln A,\ln B \right) + 
   \ln \{ 1 + \exp \left[ - |\ln A - \ln B| \right] \} $$
ought to be used.

\foilhead{\r Energy and Specific Heat Calculation}\vskip -30pt

Multicanonical data for the energy per spin (with jackknife errors) of 
the $2d$ Ising model on a $\s 20\times 20$ lattice are produced in
assignment {\tt a0501\_01} and compared with the exact results of 
Ferdinand and Fisher in assignment {\tt a0501\_03}. Results:

\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/2dI_es.pdf}
\end{center} \end{figure}

\newpage\noindent \b
The same numerical data allow to calculate the {\bf specific heat} 
defined by
$$\s  \index{multicanonical!canonical estimator!specific heat}
  C\ =\ {d\,\widehat{E}\over d\,T}\ =\ \beta^2\, \left( \,
  \langle E^2 \rangle - \langle E \rangle^2 \right)\ .$$
\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/2dI_C.pdf}
\end{center} \end{figure}

\foilhead{\r Free Energy and Entropy Calculation}\vskip -30pt
At $\s \beta=0$ the Potts partition function is $\s Z=q^N$. Therefore, 
multicanonical simulations allow for proper normalization of the 
partition function, if $\s \beta=0$ is included in the temperature 
range. Example: Entropies from multicanonical simulations of the 
The 2d Ising and 10-state Potts models on a $\s 20\times 20$ lattice 
(assignment {\tt a0501\_03} for Ising and {\tt a0501\_05} for
10-state Potts).
\begin{figure} \begin{center}
\includegraphics[height=0.58\textheight]{figures/entropy.pdf}
\end{center} \end{figure}


\foilhead{\r\Large\bf Summary}

\begin{itemize}
\item We considered Statistics, Markov Chain Monte Carlo simulations, 
      the Statistical Analysis of Markov chain data and, finally, 
      Multicanonical Sampling.

\item It is a strength of computer simulations that one can generate
      artificial (not realized by nature) ensembles, which enhance the 
      probabilities of rare events on may be interested in, or speed up 
      the dynamics.

\item Each method comes with its entire computer code.
% Much more in my book ...

\end{itemize}
\vfil
\newpage
% ------------------------------------------------------------------
%\begin{figure} \begin{center}
%\includegraphics[width=0.8\textwidth]{figures/MC3.jpg}
% \caption{\small Book cover.}
%\end{center} \end{figure}
% ------------------------------------------------------------------

\end{document}
