Let X = (x1, x2, ... , xd) be a point in a d-dimensional
space,

and dX = dx1 dx2 ... dxd the volume element.

We want to evaluate the integral over the domain Omega:

G = intwhere p(X) is interpreted as probability (density) satisfyingg(X) p(X) dX_{Omega }

p(X) > 0,The basic Monte Carlo method of calculating the integral is to draw a set (sequence) of random valuesint

p(X) dX = 1_{Omega}

X1, X2, ..., XNdistributed according to the probability p(X), and to form the arithmetic mean

GThe quantity GN is approximately G, and we say G_{N}= sumg(Xi) / N_{i=1 to N}

It is crucial to understand why p(Xi) did not multiply
g(Xi). If Xi were equally spaced in the domain Omega, we would include
p(Xi). Since Xi is distributed according to p(Xi)---**there are more points
where p(Xi) is large**---**we already take it
into account**.

*How to divide g(X) and p(X)?*

Very often the probability p(X) itself is a uniform distribution
p(X) = constant in the domain Omega. Other times, we need other form of
distribution, e.g., in statistical mechanics, we like to have Boltzmann
distribution, p(X) ~ e^{-E/kBT}. Splitting the integrand
g(X) p(X) into a probability p(X) and a quantity g(X) is somewhat arbitrary
and is for the convenience of computation. However, the choice of p(X)
does influence the statistical accuracy of the Monte Carlo estimates.

*Fundamental theorems*

We'll just state two important theorems which form the theoretical basis of the Monte Carlo method.

*(Strong) Law of Large Numbers:*

The arithmetic mean G_{N}converges with probability one to the expectation value G. In mathematical notation

P{ limThis theorem guarantees the convergence of Monte Carlo estimates to the correct answer in the limit of infinite many number of points. But it does not tell us how fast it converges. The next theorem tells us more.G_{ N to infinity}_{N}= G } = 1.

Let's define the variance as

sigmaand a normalized value (also a random variable)^{2}= int(g(X) - G)_{ Omega}^{2}p(X) dX

= < g2 > - < g >2

tthen_{N}= sqrt{N} ( G_{N}- G ) / sigma

limIn other words, GP{ a < t_{N to infinity}_{N}< b } =int

1/sqrt(2pi) e_{a to b}^{-1/2 t2}dt._{N}is a random variable distributed according to Gaussian distribution with mean G and variance sigma^{2}/N for sufficiently large N.

As N goes to infinity, the observed G_{N}turns up in ever narrower interval near G and one can predict the probability of deviations measured in units of sigma. The observed G_{N}is within one standard error (i.e. sigma/sqrt{N}) of G 68.3% of the time, within two standard errors 95.4% of the time, and within three standard errors 99.9% of the time. This is due to the property of the Gaussian distribution.

*Error of Monte Carlo calculation*

From the central limit theorem, we have

G = GThe error of Monte Carlo estimates for G is_{N}+ error

| error | = epsilon ~ sigma / sqrt{N}From the definition of sigma, we can also estimate the error epsilon as well as the value G itself:

sigmaHowever, it is better to use the following unbiased estimator for sigma2:^{2}= int(g(X) - G)_{ Omega}^{2}p(X) dX

= < g2 > - < g >2

sigma2 ~Monte Carlo error decreases as 1/sqrt{N}, as the number of samples (points) N increases. Since N is proportional to the computer CPU time T, error epsilon is proportional to TN/(N - 1){1/N sum

g_{i=1 to N}(Xi) - (1/N sum^{2}g(Xi))_{i=1 to N}}^{2}

Let's compare the efficiency of Monte Carlo method with
that of deterministic method of quadrature (e.g. trapezoidal or Simpson
rule). The deterministic methods typically have error

epsilon ~ hdue to discretization, where h is grid size. The computer time goes as^{k}

T ~ hin d dimensions when the grid size is decreased. Thus error goes as^{-d}

epsilon ~ 1/TFor Simpson rule of numerical integration, k = 4. We see that for a one-dimensional integral, the error decreases much faster than that of Monte Carlo method. However as the dimensions of the integral increase, the Simpson rule becomes worse than the Monte Carlo method for a fixed amount of computer time.^{k/d}

In conclusion, standard numerical quadrature is very good
for low dimensions. Monte Carlo method is superior for higher dimensional
integrations.

*Some examples*

Example 1:

Evaluate the integralExample 2:

intSolution:sqrt{1 - x_{0 to 1}^{2}} dx = pi/4.Let's take the probability distribution as

1 , 0 < x < 1so that x is uniformly distributed in the unit interval. And

p(x) = {

0 , otherwise

g(x) = sqrt{1 - x2}we can cast the integral in the required form

int g(x) p(x) dx.The Monte Carlo estimator for the integral is

GThis can be implemented in C with a few lines:_{N}= 1/N sumg(xi)_{i=1 to N}

= 1/N sumsqrt{1 -xi_{i=1 to N}^{2}}.

sum = 0;

for(i = 0; i < N; ++i)

sum += sqrt(1 - pow(drand48( ), 2));

mean = sum / N;

Evaluate the integral

intand estimate the Monte Carlo error.e_{0 to infinity}^{-x}cos x dx

Solution:

Since the integral domain is from 0 to infinity, it is impractical to sample the value x uniformly from 0 to infinity. However, it is quite natural to sample x from 0 to infinity by the exponential distribution:

p(x) = eand then^{-x}

g(x) = cos(x)Even though x can take any positive value, but the probability that it takes large value is very small. We sample the exponential distribution by

x = - ln xiwhere xi is a random number uniformly distributed between 0 and 1.A Monte Carlo estimator for the integral is then

GQuestion: why do we omit the exponential factor?_{N}= 1/N sumcos( ln xi_{i=1 to N}_{i})We can also estimate the error sigma/sqrt{N} in G

_{N}by a Monte Carlo calculation of sigma:

sigma^{2}~N/(N - 1) {1/N sum

cos2( ln xi_{ i=1 to N}_{i }) - G_{N}^{2}}

Homework