10.7 Monte-Carlo Integration
165
Fig. 10.2 Monte-Carlo approximation to the integral a function of number of random numbers
used
The first line generates N random numbers between −1 and 1 and second line
sums up all the function evaluations. The disadvantage is the slow convergence,
especially when integrating functions that are small over large sub-regions of the
integration interval. Integrals of Gaussians over infinite intervals are an example that
is problematic. We therefore need to find a better method; one that pays more attention
to the regions where the function is large, or even better, generate random numbers
whose distribution mimics the integrand. In other words, we try to build a random
number generator that produces random numbers, whose histogram reproduces the
function f (x). This procedure of building a tailor-made random number generator
is called importance sampling.
One way of generating random numbers that mimic the integrand is the
acceptance-rejection method which is based on generating two random numbers
x j and y j , where y j must lie between the minimum and maximum of the function
f (x). We then select random numbers by taking those x j for which y j < f (x j ). A
visualization of this method is based on throwing darts onto a target plane with x
and y and the line f (x) drawn on the target. If we hit below the line we pick the
value of x, if it is above the line, we ignore it. If we repeat this procedure sufficiently
often, the random numbers x j are clustered around values where f (x) is large. The
MATLAB code to produce such a distribution is also rather compact
x=-1+2*rand(1,N); % random numbers between -1 and 1
y=rand(1,N);
% random numbers between 0 and 1
x=x(y
% select only values under f(x)
hist(x,30);
% display histogram
where the last line only serves to display the distribution of random numbers.
Figure 10.3 shows the output of running the code for N = 10
4 iterations.
But this way of generating the distribution still suffers from a large rejection rate
of random numbers and is still very inefficient in regions where f (x) is small. We
165
Fig. 10.2 Monte-Carlo approximation to the integral a function of number of random numbers
used
The first line generates N random numbers between −1 and 1 and second line
sums up all the function evaluations. The disadvantage is the slow convergence,
especially when integrating functions that are small over large sub-regions of the
integration interval. Integrals of Gaussians over infinite intervals are an example that
is problematic. We therefore need to find a better method; one that pays more attention
to the regions where the function is large, or even better, generate random numbers
whose distribution mimics the integrand. In other words, we try to build a random
number generator that produces random numbers, whose histogram reproduces the
function f (x). This procedure of building a tailor-made random number generator
is called importance sampling.
One way of generating random numbers that mimic the integrand is the
acceptance-rejection method which is based on generating two random numbers
x j and y j , where y j must lie between the minimum and maximum of the function
f (x). We then select random numbers by taking those x j for which y j < f (x j ). A
visualization of this method is based on throwing darts onto a target plane with x
and y and the line f (x) drawn on the target. If we hit below the line we pick the
value of x, if it is above the line, we ignore it. If we repeat this procedure sufficiently
often, the random numbers x j are clustered around values where f (x) is large. The
MATLAB code to produce such a distribution is also rather compact
x=-1+2*rand(1,N); % random numbers between -1 and 1
y=rand(1,N);
% random numbers between 0 and 1
x=x(y
hist(x,30);
% display histogram
where the last line only serves to display the distribution of random numbers.
Figure 10.3 shows the output of running the code for N = 10
4 iterations.
But this way of generating the distribution still suffers from a large rejection rate
of random numbers and is still very inefficient in regions where f (x) is small. We
