Monte Carlo Integration

Demonstration Worksheet by Mike May, S.J.- maymk@slu.edu

Section 15.4 - an example

> restart:

Section 15.4 discusses the Monte Carlo method of integration.

The Monte Carlo method is a technique that can only reasonable by used with a computer.

It will not be something we spend much time on. Nevertheless, it is interesting to see a worked example.

We need to start by defining the function, and the box we use the method on.

> f := x -> sqrt(1-x^2);
lowx := 0: highx := 1:
lowy := 0: highy := 1:

> f := proc (x) options operator, arrow; sqrt(1-x^2) end;

It seems good to look at the graph of the function in the rectangle.

> plot(f(x), x=lowx..highx, y=lowy..highy, axes=BOXED);

Next we need a way to pick a random point. The Maple command rand() gives a random 12 digit positive integer. To convert to a random number between a and b we use a + (b-a)*rand()/10^12.

> randx := evalf(lowx + (highx - lowx)*rand()/10^(12));
randy := evalf(lowy + (highy - lowy)*rand()/10^(12));
evalb(randy < f(randx));

Now we want to put this in a loop, try it lots of times and see how it works. By looking at the graph, we see that the correct answer is Pi/4.

> tries := 1000:
goodtries := 0:
for i from 1 to tries do
randx := evalf(lowx + (highx - lowx)*rand()/10^(12));
randy := evalf(lowy + (highy - lowy)*rand()/10^(12));
if (randy <f(randx)) then
goodtries := goodtries+1:
fi:
od:
prob := evalf((highx - lowx)*(highy - lowy)*goodtries/tries);
answer := evalf(Pi/4);

To make it easier for you to modify the problem and experiment, all the needed commands are repeated in one block.

> f := x -> sqrt(1-x^2);
lowx := -1: highx := 1:
lowy := 0: highy := 1:
tries := 1000:
goodtries := 0:
for i from 1 to tries do
randx := evalf(lowx + (highx - lowx)*rand()/10^(12));
randy := evalf(lowy + (highy - lowy)*rand()/10^(12));
if (randy <f(randx)) then
goodtries := goodtries+1:
fi:
od:
prob := evalf((highx - lowx)*(highy - lowy)*goodtries/tries);
answer := evalf(Pi/2);

>