1. Introduction
In this exercise we investigate multi-dimensional Monte Carlo
integration. A classical example of the
Monte Carlo integration approach is to determine the area of a circular
"pond" by
throwing stones into it.
If you know:
- the total area Atot
over which you are throwing the stones, and
- that the locations (including the pond and its surrounding
area) where they landed are uniformly and randomly distributed
then the area of the pond is determined as the ratio of the number of
stones that hit the pond compared to the total number, multiplied by
the total area
Apond
= (Nhit/Ntot )* Atot
Multi-dimensional Monte Carlo integration implements this by throwing random
numbers into a multi-dimensional space and determining when the values
of these random numbers fall inside or outside of the boundary of the
volume.
For example, to determine the area of a circle (we denote this by V
(2)
, the generalized "volume" of the 2-sphere), we can
express it as the integral of an area element such as an annular ring:
or we can use a Monte Carlo approach where we use the
following
algorithm:
- generate two uniform random deviates x and y
over
an
interval [0,1]
- calculate the Cartesian distance r to the point they
correspond
to, r2 = x2 + y2
- test if r falls within the circle determined by the maximum
radius
of integration R
- if r <= R, count this test as a "hit"
- if r>R count the test as a "miss"
- After N total trials, the area of the circle is
approximated
by (Nhit / N) * 4R2
where
4R2 is the area of the box with sides
D=2R
in length which circumscribes the circle.
Simple example: area of the unit circle
Write your own version of the program (which you will use as a template
to modify later), which demonstrates how the above algorithm works on
the unit circle.
Using your favorite plotting program (identifying
those points inside versus outside by either their marker, color, or
both), produce a display like that
shown in Figure 1. Include this plot in your report, along with a
meaningful
description of what is being plotted.
Figure 1: A plot of the data generated by
the circle "volume" calculation program (plot via gnuplot).
Exercise Part I:
- Copy your program and modify it to calculate the volume of a
sphere (V(3)) using the Monte Carlo
method. Compare your result to the true volume of a sphere, quantitatively.
- Generate several different runs of your program as a function of
the number of trials N for N=10,30,100,300,1000,3000,10000,30000 and
plot the fractional error vs. number of trials.
- Make a 3-D plot of this data and include it in your report.
- Make sure to evaluate the errors and compare with statistical
expectations.
Multi-dimensional integration where you might really want Monte Carlo: The Volume of a Hypersphere
Of course the example above is kind of trivial, since we can do the
integral analytically. So we will try something
more interesting.
There are many applications in physics where integrations must be
performed
over many more than even 3-dimensions. For example,
- phase space integrals over position & momentum (dx
dy
dz dpx dpy dpz ) are
6-dimensional
- some extended quantum field theories can have 10 or more
dimensions
- Multi-dimensional integrals also occur in experiments, where the
dimensions
can include things like spatial and temporal detector response
functions
(usually without analytical forms--they must be numerically
integrated),
and very anisotropic arrival distributions for the signals of interest.
So in many cases, Monte Carlo methods are the only alternative.
If
you stay in physics (or even if not--such techniques are now widely
used
in financial modeling) you will most certainly need to become familiar
with
Monte Carlo methods.
You are now going to make a copy of your program and generalize in in
the following way.
The
problem is to calculate
the volume of an N-dimensional sphere:
a
hypersphere. We might
do this analytically if we knew the surface
area element, because we could then use an integral like the one for
the
circle above, which integrates an annular element (the boundary) of
dimension one less than
the one we want.
But what is the
3-D annular element for a 4-D sphere??
Since this is not something we know off the top of our head, we can
just
proceed to use the Monte Carlo method to compare the number of "hits"
we
get within the sphere, to the total number of trials, and then multiply
this fraction by the total volume of a 4-D
hypercube:
V(4)hypercube = D4
which we can easily
calculate even if it is hard for us to visualize more than a
3-dimensional cube.
It is of course also possible to make an analytical estimate
using integral
calculus,
but it is not at all trivial; see (
here
) for an example of how to do it.
Exercise Part II:
- Complete your program so that it
works for calculating volumes for any hypersphere dimension up to 10,
and compares
its answer with posted analytical values.
- To get a feel for what is going on,
calculate the volume of 2 and
3-dimensional spheres and check the accuracy vs. number of trials Nmax.
Show a plot of fractional error
vs. N for 2D and 3D.
- Using enough trials to get a few %
accuracy, calculate the volume
for
the unit N-hypersphere from 2 to 10 dimensions
using your program. Plot the value as a function of the number of
dimensions and note the
behavior. Is it surprising?
- Compare the value you get at each
dimension
(up
to D=7) to the
analytical values.
- Check the precision you get as a
function of Nmax and the number
of
dimensions, again up to 7-D. How many trials would
you
need to get a precision of 10-5 in 10-dimensions?