Monte-Carlo Simulations

Prerequisites: Statistics. Full path (2 pages)

• Intuitive
• Concrete
• Abstract
• History

Intuitive

Monte Carlo methods compute quantities by random sampling. A standard first example: draw a square, inscribe a quarter circle, and throw darts at random at the square. The fraction of darts that land inside the quarter circle approaches the ratio of the two areas, which is a quarter of pi. Counting hits gives an estimate of pi, and the estimate improves the more darts we throw.

The same idea works for any average or integral. Instead of evaluating a function on a regular grid of points, we evaluate it at randomly chosen points and take the mean. A grid becomes hopeless in many dimensions: with ten points per axis, a system with 100 variables would need ten to the power 100 grid points. The error of a random estimate shrinks with the number of samples at a rate that does not depend on the number of dimensions. This is why Monte Carlo is the standard tool for problems with very many variables, such as a magnet made of many interacting spins, a gas of many molecules, or a quantum field on a lattice.

In physics we usually need samples from a particular distribution, for example configurations weighted by the Boltzmann factor. Monte Carlo algorithms such as the Metropolis algorithm do this with a random walk through configuration space that visits each configuration with the right frequency, so that averages over the walk become thermal averages.

In Monte-Carlo simulations parameters are chosen randomly and then compared to a predefined goal.

Concrete

Monte Carlo integration. Suppose we want the expectation value of $f$ under a probability density $p$, $$I = \int f(x)\,p(x)\,dx.$$ Draw $N$ independent samples $x_1,\dots,x_N$ from $p$ and form the sample mean $$\hat I_N = \frac1N\sum_{i=1}^N f(x_i).$$ Its expectation value is $I$, and its variance is $\sigma_f^2/N$ with $\sigma_f^2=\int (f-I)^2p\,dx$. The statistical error is therefore $\sigma_f/\sqrt N$, independent of the dimension of $x$, and $\sigma_f$ can itself be estimated from the samples. In the dart example, $x$ is uniform on the unit square, $f=1$ inside the quarter circle and $0$ outside, so $I=\pi/4$ and $\sigma_f^2=\frac\pi4(1-\frac\pi4)$.

The error can be reduced by importance sampling: draw from a density $q$ that is large where $|fp|$ is large, and average $f(x)p(x)/q(x)$ instead. The mean is unchanged and the variance can be much smaller.

Markov chain Monte Carlo. In statistical mechanics we need averages in the canonical ensemble, $$\langle A\rangle = \frac1Z\sum_C A(C)\,e^{-\beta E(C)},\qquad Z=\sum_C e^{-\beta E(C)},$$ where $C$ runs over configurations, for example all $2^N$ spin configurations of an Ising model. We cannot sample $e^{-\beta E}/Z$ directly because $Z$ is unknown. The Metropolis algorithm produces a chain of configurations instead:

1. From the current configuration $C$, propose $C'$, for example by flipping one randomly chosen spin, with a symmetric proposal probability.
2. Compute $\Delta E = E(C')-E(C)$.
3. Accept $C'$ with probability $\min(1,e^{-\beta\Delta E})$; otherwise keep $C$ as the next element of the chain.

Only energy differences enter, so $Z$ is never needed. The acceptance rule satisfies detailed balance, $e^{-\beta E(C)}P(C\to C') = e^{-\beta E(C')}P(C'\to C)$, so the Boltzmann distribution is stationary. After discarding an initial equilibration period, $\langle A\rangle$ is estimated by the average of $A$ along the chain. Successive configurations are correlated, so the error is larger than for independent samples, and it has to be estimated with the autocorrelation time or by binning.

In lattice QCD and other lattice field theories, the same method samples field configurations with weight $e^{-S_E}$, where $S_E$ is the Euclidean action obtained after a Wick rotation.

• An introduction to Monte Carlo methods by Jean-Charles Walter, Gerard Barkema
• The Monte Carlo method in quantum field theory by Colin Morningstar

Abstract

Let $(\Omega,\mu)$ be a probability space and $f\in L^2(\mu)$. For independent samples $X_i\sim\mu$, the strong law of large numbers gives $\hat I_N=\frac1N\sum_{i=1}^Nf(X_i)\to\int f\,d\mu$ almost surely, and the central limit theorem gives $$\sqrt N\,\big(\hat I_N-\textstyle\int f\,d\mu\big)\;\xrightarrow{d}\;\mathcal{N}(0,\sigma_f^2).$$ The rate $N^{-1/2}$ is independent of the dimension of $\Omega$. Deterministic quadrature rules on a grid in $d$ dimensions converge as $N^{-k/d}$ for integrands of smoothness $k$, so in high dimension sampling wins. Quasi-Monte Carlo methods replace random points by low-discrepancy sequences and can reach errors close to $N^{-1}$ for sufficiently regular integrands.

Markov chain Monte Carlo constructs a transition kernel $P(x,dy)$ that leaves the target $\pi$ invariant, $\int\pi(dx)P(x,dy)=\pi(dy)$. Detailed balance, $\pi(dx)P(x,dy)=\pi(dy)P(y,dx)$, is a sufficient condition and makes $P$ self-adjoint on $L^2(\pi)$. For a Markov chain on a finite state space that is irreducible and aperiodic, the distribution of $X_n$ converges to $\pi$ from any starting point, and ergodic averages converge to $\pi$-expectations; analogous statements hold on general state spaces under Harris recurrence. The Metropolis-Hastings rule turns any proposal kernel $q(x,y)$ into such a chain by accepting with probability $$\alpha(x,y)=\min\left(1,\frac{\pi(y)\,q(y,x)}{\pi(x)\,q(x,y)}\right),$$ which requires $\pi$ only up to normalization.

For a stationary chain the variance of the ergodic average is $$\mathrm{Var}(\hat I_N)\approx\frac{2\tau_{\mathrm{int}}\,\sigma_f^2}{N},\qquad \tau_{\mathrm{int}}=\frac12+\sum_{t=1}^\infty\rho_f(t),$$ where $\rho_f(t)$ is the normalized autocorrelation of $f(X_n)$ at lag $t$, so $N/(2\tau_{\mathrm{int}})$ is the effective number of independent samples. The decay of autocorrelations is controlled by the spectral gap of $P$. Near a continuous phase transition $\tau_{\mathrm{int}}$ grows as a power of the correlation length for local update algorithms (critical slowing down), which motivates nonlocal updates such as cluster algorithms for spin models; lattice gauge theory with dynamical fermions uses hybrid Monte Carlo, which updates all variables at once along molecular dynamics trajectories. All of this assumes a nonnegative weight. When the Boltzmann weight is not positive, as for fermions at finite density or real-time path integrals, the reweighted averages suffer from cancellations that grow exponentially with the volume, the sign problem.

History

• The Beginning of the Monte Carlo Method by N. Metropolis

Path to this page

1. Calculus
2. Statistics