Integrals 215 problem 14 (Monte Carlo).mws

1. Monte Carlo integration.

The basic idea of Monte Carlo integration is to use probability to estimate the value of an integral. Suppose that you have a curve y = f(x) , which for simplicity we will assume is nonnegative, and you want to find the area underneath the curve from x = a to x = b . Choose some number M such that 0 <= f(x) <= M for all x in the interval

[a, b] . Pick n points [x[i], y[i]] at random in the rectangle bounded by [a, b] on the x -axis and [0, M] on the y -axis. Count how many of these points satisfy y[i] <= f(x[i]) . If there are m such points, then estimate that the value of the integral is m*A/n , where A = (b-a)*M is the area of the rectangle containing the curve..

Since the purpose of the activities here is to apply technology to solve problems in mathematics, you do not have to pay close attention to the manner in which the function below is defined. What is important is that you know how to use the function to do Monte Carlo integration. That we will describe below after defining a Maple proc that will do the random selection of points and the evaluation of the Monte Carlo estimate of the integral for us.

> mcint:=proc(f,a,b,M,n)

> local i,x,y, MM, hits, pointsin, rf, pointsout,pin,pout,pf,actual,estimate,err,percenterror;
hits:=0;

> MM:=100000*n;

> pointsin:=[];

> pointsout:=[];

> rf:=rand(0..MM);

> for i to n do

> x:=evalf(a+(b-a)*rf()/MM);

> y:=evalf(M*rf()/MM);

> if (y<=f(x)) then

> hits:=hits +1;

> pointsin:=[op(pointsin),[x,y]];

> else

> pointsout:=[op(pointsout),[x,y]];

> fi;

> od;

> pin:=plots[listplot](pointsin,style = POINT, symbol=CIRCLE,color=blue):

> pout:=plots[listplot](pointsout,style=POINT,symbol=DIAMOND,color=red):

> pf:=plot(f('x'),'x'=a..b,'y'=0..M, color=GREEN):

> print(plots[display]({pout,pin,pf}));

> actual := evalf(int(f(z), z=a..b));

> estimate := evalf((b-a)*M*hits/n);

> err:= estimate-actual;

> percenterror := err*100/actual:
printf("Actual Value: %f\n",actual);
printf("Estimated Value: %f\n",estimate);

> printf("Relative error: %f%%\n",percenterror);

> end;

mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...
mcint := proc (f, a, b, M, n) local i, x, y, MM, hi...

Let us choose a simple function to perform Monte Carlo integration on. The example has the nice feature that we can compute the exact value of the integral by hand. Of course, in practice, we would want to use the Monte Carlo method on functions that we cannot integrate exactly. It should be clear to the reader that in the input to the mcint proc, f is the function, a and b are the left and right endpoints of the interval, M is a number which should be chosen by the reader to be greater than or equal to the maximum value of the function on the interval [a,b] and n is the number of points to use for the estimate. Because the routine will print a picture showing the points chosen, I would not recommend using too large a value for n.

Try executing the routine below several times. You will notice that the estimate of error changes each time.

> f:=x->x^2;a:=0;b:=1;M:=2;n:=1000;

f := proc (x) options operator, arrow; x^2 end proc...

a := 0

b := 1

M := 2

n := 1000

> mcint(f,a,b,M,n);

[Maple Plot]

Actual Value: .333333

Estimated Value: .330000

Relative error: -1.000000%

Submission:

For each of the integrals below, do the following.

a) Define the function and the endpoints a and b in Maple.

b) Determine a reasonable value of M.

c) Apply the Monte Carlo routine to get an estimate of the integral.

1) Int(x*exp(-2*x),x = 0 .. 1) .

>

2) Int(ln(x)^3,x = 1 .. 2) .

Submission worksheet:

>

>

>