<Text-field style="_pstyle256" layout="_pstyle256">Games of Chance: Exploring the Monte Carlo Technique for Numerical Integration and Computing Probabilities with Improper Integrals</Text-field>

<Text-field style="_pstyle256" layout="_pstyle256">Games of Chance: Exploring the Monte Carlo Technique for Numerical Integration and Computing Probabilities with Improper Integrals</Text-field>

Note: You may notice differences between this Maple worksheet and the equivalent Mathematica notebook. These differences were introduced to preserve the content of these modules and were necessary because of major functional differences between Maple and Mathematica.

<Text-field style="Heading 1" layout="Heading 1">Introduction</Text-field>

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.

<Text-field style="_cstyle259" layout="Heading 1"><Font size="18">Technology Guidelines</Font></Text-field>

NOTE: If you have just finished a worksheet, restart Maple before executing a new worksheet. TO OPEN SECTIONS, Click on the Arrow sign at the left hand side of the screen or select Expand All Sections from the View drop down menu.

TO STOP AN EXECUTION Click on STOP button from the toolbar.

ORDER OF EXECUTION Execute commands in the order given. Do not skip any Maple Input lines within a given worksheet

Alternatively, you can execute the entire worksheet by selecting the Execute Worksheet command from the Edit drop down menu.

SAVING WORKSHEETS.

You can save anytime to any directory you choose, and it is wise to save often. EXPERIENCING MAJOR PROBLEMS Save if appropriate, and then shut down Maple and start it up again.

<Text-field style="Heading 1" layout="Heading 1">Part I: The Monte Carlo Method</Text-field>

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 used here as an example cannot be integrated symbolically, so we are using Maple'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.

restart;

f := x->(1/(sqrt(2*Pi))*exp(-0.5*x^2));

xmin := -2:

xmax := 2:

ymin := 0:

ymax := 1:

pf := plot(f(x), x=xmin..xmax, y=ymin..ymax, color=GREEN):

pf;

In utilizing the Random command in Maple, 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.

numberofpoints := 500:

areaofbox := (xmax-xmin)*(ymax-ymin):

count := 0:

listout := []:

listin := []:

rf := rand(0..10000):

for i from 1 to numberofpoints do px := evalf(xmin + (xmax-xmin)*rf()/10000): py := evalf(ymin + (ymax-ymin)*rf()/10000): if py < evalf(f(px)) then count := count+1: listin := [op(listin), [px,py]]: else listout := [op(listout), [px,py]]: fi: od:

pout := plots[listplot](listout,style = POINT, symbol=CIRCLE):

pin := plots[listplot](listin, style=POINT, symbol=POINT, color=RED):

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

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

plots[display]({pout,pin,pf}); actual := evalf(int(f(x), x=xmin..xmax)); estimate := evalf(areaofbox*count/numberofpoints); err:= estimate-actual; percenterror := err/actual: print(`relative error = `,percenterror);

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? Analyze the error involved in such an approximation. In a later module on using Monte Carlo techniques for approximating volumes, 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.

<Text-field style="Heading 1" layout="Heading 1">You Try It: Part I</Text-field>

<Text-field style="Heading 1" layout="Heading 1">Generating 1000 Random Points</Text-field>

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

numberofpoints := 500:

areaofbox := (xmax-xmin)*(ymax-ymin):

count := 0:

listout := []:

listin := []:

rf := rand(0..10000):

for i from 1 to numberofpoints do px := evalf(xmin + (xmax-xmin)*rf()/10000): py := evalf(ymin + (ymax-ymin)*rf()/10000): if py < evalf(f(px)) then count := count+1: listin := [op(listin), [px,py]]: else listout := [op(listout), [px,py]]: fi: od:

pout := plots[listplot](listout,style = POINT, symbol=CIRCLE):

pin := plots[listplot](listin, style=POINT, symbol=POINT, color=RED):

plots[display]({pout,pin,pf}); actual := evalf(int(f(x), x=xmin..xmax)); estimate := evalf(areaofbox*count/numberofpoints); err := estimate-actual; percenterror := err/actual: print(`relative error = `,percenterror);

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

<Text-field style="Heading 1" layout="Heading 1">Selecting a Different Function</Text-field>

When applying a function of your own to this procedure, be certain to determine the range in which you enclose your desired area (values of xmin, smax, 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 your function, minimum and maximum values for your variables, except for leaving your ymin at 0o. You should realize that the bigger your rectangular region, the more points you will need to generate to get an answer that is reasonably accurate. Why?

restart; f := x->1-sin(x^2);

xmin := -3:

xmax := 3:

ymin := 0:

ymax := 2:

pf := plot(f(x), 'x'=xmin..xmax, y=ymin..ymax, color=GREEN): print(pf);

numberofpoints := 500:

areaofbox := (xmax-xmin)*(ymax-ymin):

count := 0:

listout := []:

listin := []:

rf := rand(0..10000):

for i from 1 to numberofpoints do px := evalf(xmin + (xmax-xmin)*rf()/10000): py := evalf(ymin + (ymax-ymin)*rf()/10000): if py < evalf(f(px)) then count := count+1: listin := [op(listin), [px,py]]: else listout := [op(listout), [px,py]]: fi: od:

pout := plots[listplot](listout,style = POINT, symbol=CIRCLE):

pin := plots[listplot](listin, style=POINT, symbol=POINT, color=RED): plots[display]({pout,pin});

plots[display]({pout,pin,pf}); actual := evalf(int(f(x), x=xmin..xmax)); estimate := evalf(areaofbox*count/numberofpoints); err := estimate-actual; percenterror := err/actual: print(`relative error = `,percenterror);

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?

<Text-field style="Heading 1" layout="Heading 1">Part II: Probability Density Functions and Improper Integrals</Text-field>

<Text-field style="_cstyle256" layout="Heading 1"><Font size="18">The Standard Normal Distribution</Font></Text-field>

The function that we looked at above LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2J0YrLUYjNiYtRiw2JFEiZkYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y3USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQS8lKXN0cmV0Y2h5R0ZBLyUqc3ltbWV0cmljR0ZBLyUobGFyZ2VvcEdGQS8lLm1vdmFibGVsaW1pdHNHRkEvJSdhY2NlbnRHRkEvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZQLUkobWZlbmNlZEdGJDYlLUYjNiQtRiw2JFEieEYnRjZGNi8lJ2l0YWxpY0dGQUY9RjYtRjo2LlEiPUYnRmVuRj1GP0ZCRkRGRkZIRkpGTC9GT1EsMC4yNzc3Nzc4ZW1GJy9GUkZbby1JJm1mcmFjR0YkNigtRiM2Ji1JI21uR0YkNiVRIjFGJ0ZlbkY9LUY6Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y9Rj9GQkZERkZGSEZKRkxGTkZRLUklbXN1cEdGJDYlLUY6Ni5RLyZFeHBvbmVudGlhbEU7RidGZW5GPUY/RkJGREZGRkhGSkZMRk4vRlJRLDAuMTExMTExMWVtRictRiM2Ji1GOjYuUSomdW1pbnVzMDtGJ0ZlbkY9Rj9GQkZERkZGSEZKRkwvRk9RLDAuMjIyMjIyMmVtRicvRlJGZ3AtRiM2Jy1GY282JVEkMC41RidGZW5GPUZmby1GIzYkLUZqbzYlRlgtRmNvNiVRIjJGJ0ZlbkY9LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y2RitGNkYrRjZGZXFGNi1GIzYkLUkmbXNxcnRHRiQ2Iy1GIzYmRmJxRmZvLUYsNiVRJSZwaTtGJ0ZlbkY9RjZGNi8lLmxpbmV0aGlja25lc3NHRmVvLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmZyLyUpYmV2ZWxsZWRHRkFGNkYrRjY= , 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 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JS1JI21vR0YkNi1RKiZ1bWludXMwO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGOi8lKXN0cmV0Y2h5R0Y6LyUqc3ltbWV0cmljR0Y6LyUobGFyZ2VvcEdGOi8lLm1vdmFibGVsaW1pdHNHRjovJSdhY2NlbnRHRjovJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0ZJLUYsNiVRKCZpbmZpbjtGJy8lJ2l0YWxpY0dRJXRydWVGJy9GNlEnaXRhbGljRidGNUYrRjU= to +LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJmluZmluO0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=, 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.

restart; f := x->1/(sqrt(2*Pi))*exp(-0.5*x^2);

plot(f(x), x=-3..3); print(`total mass = `, int(f(x), x=-infinity..infinity)); print(`mean = `, int(x*f(x), x=-infinity..infinity)); var := int(x^2*f(x), x=-infinity..infinity): print(`variance = `, var); print(`standard deviation = `, sqrt(var));

<Text-field style="Heading 1" layout="Heading 1">Computing Probabilities for the Normal Distribution</Text-field>

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

print(`probability that x is greater than .5 is `, int(f(x), x=.5..infinity));

Similarly, to find the probability that x is less than -.67, we find the area under the density function from -LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJmluZmluO0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= to -.67 .

print(`probability that x is less than -0.67 is `, int(f(x), x=-infinity..-0.67));

Therefore the "-.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?

<Text-field style="Heading 1" layout="Heading 1">You Try It: Part II</Text-field>

What if you wished to examie 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 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2KEYrLUYjNiYtRiw2JFEiZkYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y3USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQS8lKXN0cmV0Y2h5R0ZBLyUqc3ltbWV0cmljR0ZBLyUobGFyZ2VvcEdGQS8lLm1vdmFibGVsaW1pdHNHRkEvJSdhY2NlbnRHRkEvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZQLUkobWZlbmNlZEdGJDYlLUYjNiQtRiw2JFEidEYnRjZGNi8lJ2l0YWxpY0dGQUY9RjYtRjo2LlEiPUYnRmVuRj1GP0ZCRkRGRkZIRkpGTC9GT1EsMC4yNzc3Nzc4ZW1GJy9GUkZbby1GIzYmLUYsNiVRKSZsYW1iZGE7RidGZW5GPS1GOjYtUTEmSW52aXNpYmxlVGltZXM7RidGPUY/RkJGREZGRkhGSkZMRk5GUS1JJW1zdXBHRiQ2JS1GOjYuUS8mRXhwb25lbnRpYWxFO0YnRmVuRj1GP0ZCRkRGRkZIRkpGTEZOL0ZSUSwwLjExMTExMTFlbUYnLUYjNiYtRjo2LlEqJnVtaW51czA7RidGZW5GPUY/RkJGREZGRkhGSkZML0ZPUSwwLjIyMjIyMjJlbUYnL0ZSRmNwLUYjNiZGX29GYm9GWEY2RitGNi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGNkYrRjZGK0Y2, for certain values of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEpJmxhbWJkYTtGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnRjI=. 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.

restart; lambda := .25: f := t->lambda*exp(-lambda*t): plot(f(x),x=0..10, labels=["time","probability"]); print(`total mass = `, int(f(x), x=0..infinity)); mean := int(x*f(x), x=0..infinity): print(`expected wait = `, mean); var := int((x-mean)^2*f(x), x=0..infinity): print(`variance = `, var); 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 less than 3 minutes Change the numbers to the appropriate values for each of the cases above.

int(f(x), x=1..5);

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=

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