Monte Carlo Technique for Numerical Integration

Adapted from material in Thomas' Calculus (10th Edition) by Finney, Weir, and Giordano

Email completed labs to tleise@amherst.edu by noon next Wednesday.

Introduction

How can you use a game of chance to evaluate a multiple integral? That is just what you will do with this lab. You will generate random points within a fixed region and then estimate the area or volume of the desired portion by considering the proportion of random points that fall within the boundaries of the desired portion. Since this will only estimate the integral, you will also explore the accuracy of this method.

Consider a game where you throw darts at a one foot by one foot square board with a circle of radius 6 inches drawn on it.  If you have absolutely no skill at darts and hit points on the board randomly, then you would expect the proportion of darts landing inside the circle would be equal to the area inside the circle (given that the area of the square is 1).  The code below throws n darts at the board randomly, counts how many hit, then plots the dart board with the hits on it.  Random[Real,{0,1}] randomly picks a real number between 0 and 1.

board=Graphics[{Hue[0],Rectangle[{0,0},{1,1}],Hue[0.6],Disk[{1/2,1/2},1/2]}];  (* EVALUATE THIS CELL -- HIT SHIFT-ENTER *)

n=100;
count=0;
darts={};
Do[{Clear[x,y],x=Random[Real,{0,1}],y=Random[Real,{0,1}],
		      If[(x-1/2)^2+(y-1/2)^2<1/4,count=count+1],AppendTo[darts,{x,y}]},{n}]
dartsplot=ListPlot[darts,PlotStyle->{PointSize[0.015],Hue[.2]}];
Show[board,dartsplot,AspectRatio->Automatic]
Print[count," darts landed inside the circle"]

Exercise 1. Estimate the area of the circle using this "darts" idea.  For example, if 82 out of 100 darts land inside the circle, we would estimate the area of the circle to be 0.82.  Get three different estimates for n=1000 darts.  How do they compare to the actual area?  Calculate the percent error = 100×(estimate-actual)/actual for each estimate.

The Monte Carlo Method

The Monte Carlo method is an example of a numerical integration technique. In this method, the volume under a surface is completely enclosed in a rectangular box and points within that box are randomly selected. The number of points that land under the surface is determined and the volume under the surface to the bottom of the box is approximated as the proportion of the generated points lying under the surface times the volume of the entire box.  Note that the sample function used here cannot be integrated symbolically, so we are using Mathematica's numerical integrator, which uses a different approximation technique, to come up with an answer that we would expect to be close to the actual one.

Our goal is to estimate the volume under the surface z=(e^(-2(x^2)-(y^2)))for 0≤x≤2, 0≤y≤1.  The height of the surface is always less than 1 and greater than 0, so the region under the surface lies inside the box [0,2]×[0,1]×[0,1].

f[x_,y_]:=Exp[-2(x^2)-(y^2)]; (* defines a function of x and y *)
surface=Plot3D[f[x,y],{x,0,2},{y,0,1},PlotRange->All]            (* EVALUATE THIS CELL -- HIT SHIFT-ENTER *)
actual=NIntegrate[f[x,y],{x,0,2},{y,0,1}];
Print["Actual volume is approximately ",actual]

Now we will generate random points in our specified region and check to see whether or not they are under our surface. The points in listin are under the surface in the volume we are computing and they will be plotted in red. The points in listout are not under the surface and they will be plotted in black. The ScatterPlot command requires the Graphics3D package which must be read in first.

n=100;
volumeofbox=2;
count=0;
listout={};
listin={};
Do[{Clear[x,y,z],x=Random[Real,{0,2}],y=Random[Real,{0,1}],z=Random[Real,{0,1}],
		      If[z<f[x,y],{count=count+1,AppendTo[listin,{x,y,z}]},AppendTo[listout,{x,y,z}]]},{n}]     (* EVALUATE THIS CELL -- HIT SHIFT-ENTER *)
plotout=ListPointPlot3D[listout];
plotin=ListPointPlot3D[listin,PlotStyle->RGBColor[1,0,0]];
Show[plotin,plotout,surface,PlotRange->All]
Print[count," points were under the surface out of ",n," total points"]

Exercise 2. Estimate the volume under the surface z=f(x,y) by using n=1000 points. Remember to take into account the volume of the total box.  How does this estimate compare to the actual volume?  Calculate the percent error of this estimate.

Error Analysis

Each time we do a Monte Carlo calculation for an integral, we can get a slightly different answer.  Now we want to repeat the Monte Carlo procedure many times to look at the pattern of the errors that result.  This gives us an idea of the "typical" error to expect as well as the range of possible errors that are likely to occur.

Clear[x,y,z]
n=1000;
errorlist={};  
repeat=200;
count=0;
Do[{count=0,Do[{Clear[x,y,z],x=Random[Real,{0,2}],y=Random[Real,{0,1}],z=Random[Real,{0,1}],
		      If[z<f[x,y],count=count+1]},{n}],
		      estimate=N[volumeofbox  count/n],error=100(estimate-actual)/actual,AppendTo[errorlist,error]},
	{repeat}]
ListPlot[errorlist, Joined->True,PlotLabel->"Percentage Errors for 200 Trials",AspectRatio->1/6]

These are referred to as "random errors," and they should be centered around 0 and should show no particular pattern. We can order them to see the range of errors that occur, as well as a histogram to see how often different sizes of errors occur.

sortlist=Sort[errorlist]
ListPlot[sortlist, Joined->True]
Histogram[errorlist]

Exercise 3. Repeat this error analysis procedure, changing the number of points, say, to n=500, n=5000, and then n=20000, to see what happens to the errors (leave the number of trials at 200, “repeat=200”).  Be patient with n=20000, it may take few minutes.  Please type a sentence summarizing the results.