Games of Chance:
Exploring the Monte Carlo Technique for Numerical Integration and Computing Probabilities with Improper Integrals

Introduction

OBJECTIVE: Learn the Monte Carlo method for approximating an integral that cannot be integrated symbolically, and see an application of improper integrals.

How can you use a game of chance to evaluate an integral, and what do games of chance have to do with improper integrals? This module will give you insight into both of these issues.
1.  You will generate random points within a fixed region and then estimate the area of the desired portion by considering the percentage of random points that fall within the boundaries of the desired portion. 
2.  You will consider a probability density function defined over an infinite interval and explore the computation of means, moments, and probabilities associated with the density function.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: The Monte Carlo Method

Following is an example of a numerical integration technique called the Monte Carlo method. In this method, the area between a curve (with nonnegative y-values) and the x-axis is completely enclosed in a rectangle, and points are randomly selected within that rectangular region. Then the number of points that land under the curve is determined, and the area under the curve (down to the x-axis) is approximated as the percentage of generated points that lie under the curve times the area of the entire rectangle. The function we use here as an example cannot be integrated symbolically, so we are using Mathematica's numerical integrator to come up with an answer that we would expect to be very close to the actual one.

The following commands will make the dimensions of the box selected 4 by 1, giving an area of 2. First we will graph our function.

Clear[x,y,f]
f[x_]:=((E^(-0.5` (x^2)))/(Sqrt[2 π]))
xmin=-2;
xmax=2;
ymin=0;
ymax=0.5`;
pf=Plot[f[x],{x,xmin,xmax},PlotRange->{{xmin,xmax},{ymin,ymax}},AxesLabel->{x,y}]

In utilizing the Random command in Mathematica, we specify the type of numbers we want and the interval over which we want those numbers randomly selected. Next we will generate 500 points and see how they fall.

\[MathematicaIcon] About Mathematica

numberofpoints=500;
areaofbox=(xmax-xmin) (ymax-ymin);
count=0;
listout={};
listin={};
Do[{Clear[x,y],x=RandomReal[{xmin,xmax}],y=RandomReal[{ymin,ymax}],If[y<f[x],{count=count+1,AppendTo[listin,{x,y}]},AppendTo[listout,{x,y}]]},{numberofpoints}]
pout=ListPlot[listout,DisplayFunction->Identity,PlotStyle->PointSize[0.01`]];
pin=ListPlot[listin,PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity,PlotStyle->PointSize[0.03`]];
Show[pin,pout,DisplayFunction->$DisplayFunction,AxesLabel->{"x","y"}]

The red points are those that are under our curve. Do they look right?

Show[pin,pout,pf,DisplayFunction->$DisplayFunction,AxesLabel->{"x","y"}]
actual=NIntegrate[f[x],{x,xmin,xmax}];
Print["actual = ",actual]
estimate=N[((areaofbox count)/numberofpoints)];
Print["red estimate = ",estimate]
error=estimate-actual;
Print["error = ",error]
relativeerror=(error/actual);
Print["relative error = ",relativeerror]

We have given the actual value as determined by the software and have computed the error. How does your error compare to those of your classmates? Each time you run this, you will get a different answer. It might be better or it might be worse. Does this method seem very accurate?  In a later module on using Monte Carlo techniques for approximating volumes, we analyze the error involved in such an approximation . In the You Try It exercises below, see if your error gets smaller when you select more random points. You will also check out the technique for a function of your own choosing.

You Try It: Part I

Generating 1000 Random Points

Redo the Monte Carlo method above, but generate twice as many random points. Compare your error to the results above. Change the number in red.

numberofpoints=500;areaofbox=(xmax-xmin) (ymax-ymin);count=0;listout={};listin={};Do[{Clear[x,y],x=RandomReal[{xmin,xmax}],y=RandomReal[{ymin,ymax}],If[y<f[x],{count=count+1,AppendTo[listin,{x,y}]},AppendTo[listout,{x,y}]]},{numberofpoints}]
pout=ListPlot[listout,DisplayFunction->Identity,PlotStyle->PointSize[0.01`]];
pin=ListPlot[listin,PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity,PlotStyle->PointSize[0.03`]];
Show[pin,pout,pf,DisplayFunction->$DisplayFunction,AxesLabel->{"x","y"}]
actual=NIntegrate[f[x],{x,xmin,xmax}];
Print["actual = ",actual]
estimate=N[((areaofbox count)/numberofpoints)];
Print["red estimate = ",estimate]
error=estimate-actual;
Print["error = ",error]
relativeerror=(error/actual);
Print["relative error = ",relativeerror]

How does your percentage error compare to that when only 500 points were generated?

Selecting a different function

When applying a function of your own to this procedure, be certain to determine carefully the size of the rectangular box in which you enclose your desired area (values of xmin, xmax, ymin, ymax). In choosing a function, choose one that is not negative in the specified domain, and note that you are finding the area between the surface and the x-y plane. The only commands that you need to enter in the following template are the ones in red, your function and minimum and maximum values for your variables, except for leaving your ymin at 0. You should realize that the bigger your rectangular region, the more points you will need to generate (red value in the second cell) to get an answer that is reasonably accurate. Why?

Clear[x,y,f]
f[x_]:=1-Sin[(x^2)]
xmin=-3;xmax=3;ymin=0;ymax=2;pf=Plot[f[x],{x,xmin,xmax},AxesLabel->{x,y}]

numberofpoints=500;areaofbox=(xmax-xmin) (ymax-ymin);count=0;listout={};listin={};Do[{Clear[x,y],x=RandomReal[{xmin,xmax}],y=RandomReal[{ymin,ymax}],If[y<f[x],{count=count+1,AppendTo[listin,{x,y}]},AppendTo[listout,{x,y}]]},{numberofpoints}]
pout=ListPlot[listout,DisplayFunction->Identity,PlotStyle->PointSize[0.01`]];
pin=ListPlot[listin,PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity,PlotStyle->PointSize[0.03`]];
Show[pin,pout,DisplayFunction->$DisplayFunction,AxesLabel->{"x","y"}]

Show[pin,pout,pf,DisplayFunction->$DisplayFunction,AxesLabel->{"x","y"}]
actual=NIntegrate[f[x],{x,xmin,xmax}];
Print["actual = ",actual]
estimate=N[((areaofbox count)/numberofpoints)];
Print["red estimate = ",estimate]
error=estimate-actual;
Print["error = ",error]
relativeerror=(error/actual);
Print["relative error = ",relativeerror]

How does your percentage error compare with the others you've done? If it is larger, is your rectangle larger? If it is smaller, is your rectangle smaller?

Part II: Probability Density Functions and Improper Integrals

The Standard Normal Distribution

The function that we looked at above, f(x) =(1/(Sqrt[2π])) (e^(-((x^2)/2))), happens to represent the probability density function for the standard normal distribution, a distribution that is used extensively in experimentation and in probability and statistics. We will use our computational tool to find:
		the total mass, which is the area under this curve from -∞ to +∞,
		the mean value of x, which is the first moment about 0,
		the variance of x, which is the second moment about the mean, and
		the standard deviation, which is the square root of the variance.

Clear[x,y,f]
f[x_]:=((E^(-((x^2)/2)))/(Sqrt[2 π]))
pf=Plot[f[x],{x,-3,3},AxesLabel->{x,y}]
Print["total mass = ",(∫_(-∞)^(+∞))f[x]\[DifferentialD]x]
Print["mean = ",(∫_(-∞)^(+∞))x f[x]\[DifferentialD]x]
Print["variance = ",var=(∫_(-∞)^(+∞))(x^2) f[x]\[DifferentialD]x]
Print["standard deviation = ",(Sqrt[var])]

Computing Probabilities for the Normal Distribution

The bell-shaped curve above is used extensively in applied problems, including error analysis. To compute the probability that x is greater than 0.5, we would find the area under the density function from 0.5 to ∞.

Print["The probability that x is greater than 0.5 is ",prob1=Integrate[f[x],{x,.5,+∞}]]

Similarly, to find the probability that x  is less than -0.67, we find the area under the density function from -∞ to -0.67 .

Print["The probability that x is less than -0.67 is ",prob2=Integrate[f[x],{x,-∞,-.67}]]

Therefore the -0.67 value of x represents approximately the 25th percentile for this distribution. What does the term percentile mean to you? Does it correspond to this interpretation?

You Try It: Part II

What if you wished to examine the length of time that people spend waiting in line at a drive-up fastfood restaurant. Frequently, the waiting times follow what is called an exponential distribution, whose density function is written as f(t)=λ (e^(-λt)), for certain values of λ. You should note that this function will be defined only for values of t greater than 0, since it would make no sense to talk about waiting a negative amount of time. Assume that time is measured in minutes.

Clear[t,λ,f,x]
λ=0.25
f[t_]:=λ (E^(-λ t))
pf=Plot[f[t],{t,0,10},{AxesLabel->{time,probability},AxesOrigin->{0,0}},PlotRange->{{0,10},{0,f[0]}}]
Print["total mass = ",(∫_0^(+∞))f[t]\[DifferentialD]t]
Print["expected wait = ",mean=(∫_0^(+∞))t f[t]\[DifferentialD]t," minutes."]
Print["variance = ",var=(∫_0^(+∞))(((t-mean))^2) f[t]\[DifferentialD]t]
Print["standard deviation of wait time = ",(Sqrt[var])," minutes"]

Consider the shape of the exponential distribution. Observe that the horizontal axis represents waiting times, while the vertical axis represents the relative frequency of those waiting times. This distribution applies to many waiting time scenarios for which most people wait a shorter period of time and only a few wait a longer period of time.

Compute the probability of waiting more than 6 minutes.

Compute the probability of waiting less than 3 minutes.

Change the numbers in red to the appropriate values to answer each of the questions posed.

Integrate[f[x],{x,1,5}]

Now go back to the set of commands and change the value of λ (lambda) (in red). See what happens to the mean and standard deviation. What does the value of lambda mean in the context of the problem?

\[MathematicaIcon] About Mathematica

In Mathematica, plots can be suppressed when first defined and then shown later with the Show command. To do this, call on the option DisplayFunction->Identity within the command that defines the plot, and then call on the option DisplayFunction->$DisplayFunction when you are ready to Show the plot. 
Go back.