Monte Carlo Integration with TI-Nspire

Forest W. Arnold

October 2020

Typeset in L

This work is licensed under the Creative Commons Attribution-Noncommercial-Share

Alike 4.0 International License. To view a copy of this license, visit http://creativecommons.

org/licenses/by-nc-sa/4.0 or send a letter to Creative Commons, 171 Second

Street, Suite 300, San Francisco, California, 94105, USA. You can use, print, dupli-

cate, share this work as much as you want. You can base your own work on it and

reuse parts if you keep the license the same.

Trademarks

TI-Nspire is a registered trademark of Texas Instruments, Inc.

1 Introduction

Monte Carlo Methods are procedures for solving deterministic problems using proba-

bilistic methods. These methods were developed to study neutron-diffusion problems

that were impossible to solve analytically [Farlow]. Monte Carlo Integration is a pro-

cedure for approximating an integral using probability. It is used to ﬁnd approximate

solutions to higher-dimensional integrals or when traditional numerical methods for

approximating integrals are impractical or unstable.

This article describes two simple Monte Carlo methods (the ”hit or miss” technique

and the sample mean technique) for ﬁnding an approximate value for an integral and

presents examples of the techniques. The TI-Nspire documents that accompany this

article contain several user-deﬁned functions that implement the techniques.

2 The ”Hit or Miss” Monte Carlo Integration Method

Recall that the value of a deﬁnite integral equals the area between the graph of the

integrand and the x-axis:

Area =

f (x)dx

Figure 1 illustrates the relation between the value of the integral of f (x) = x

over the

interval from x = 1 to x = 3 and the area between the graph of f (x) and the x-axis.

Since the area under the graph of a function and the value of a deﬁnite integral are

Figure 1: The Value of an Integral as an Area

equal, a technique for determining the area yields the value of the integral. A (very)

simple example of such a technique is ﬁnding the integral of the function f (x) = x − 2

from x = 2 to x = 5.

First, the maximum value of the function in the interval is determined, which is y =

f (5) = 5 − 2 = 3.

Next, a rectangle enclosing the graph of the function over the interval of integration

and the maximum value of the function in the interval is constructed.

The ﬁnal step is determining the percentage of the area of the rectangle between the

graph of the function and the base of the rectangle (the x axis). Multiplying the area of

the rectangle by the percentage results in an estimate of the integral.

For this simple function, the result is a rectangle containing a right triangle as shown

in Figure 2. The graph of the function divides the rectangle into two equal triangular

Figure 2: An Integral as a Percentage of the Area of a Rectangle

areas, so the area under the graph is 50 percent of the area of the containing rectangle:

area under the graph = .5 × area o f the rectangle = .5 × 9 = 4.5

which is the same as the value of the deﬁnite integral of f (x) from x = 2 to x = 5.

The example above provides the basic approach to estimating the value of a deﬁnite

integral as the area under a graph. However, it does not say anything about how to ﬁnd

the area under an arbitrary curve. This is where probability enters.

Suppose that darts are thrown at a rectangle containing the graph of a function. Some

of the darts will hit inside the area under the graph and some will hit in the area above

the graph. Count the number that hit inside the area under the graph, then calculate

the percentage that hit the area under the graph vs the total number of darts thrown.

Multiplying this percentage times the area of the rectangle yields an estimate of the

area under the graph and thus an estimate of the value of the integral.

Throwing darts at a rectangular area can be simulated using random numbers from

a uniform distribution

to generate the random x and y values for a point inside the

rectangle. The x, y values for points inside the rectangle with width w and height h are

generated as follows:

Given r1 and r2, uniformly distributed random numbers between 0 and 1, and a equal

to the left-hand x coordinate of the rectangle, x = a + r1 × w and y = r2 × h.

2.1 ”Hit or Miss” Monte Carlo Integration with TI-Nspire

Key observations for performing ”hit or miss” integration with TI-Nspire are:

1. The TI-Nspire function rand() generates uniformly distributed random num-

bers between 0 and 1.

2. The width of the rectangle is b − a, where a,b are the lower and upper bounds of

the integral.

3. The height of the rectangle is the maximum value of the integrand in the interval

[a, b]. For f (x) and f (x, y), this can be determined graphically or using calcu-

lus. For functions with more than 2 variables, calculus techniques are generally

required.

4. To determine if the random point (rx,ry) is under the graph of the integrand f (x),

compare ry and f (rx). If ry < f (rx), then the point (rx, ry) is under the graph of

the integrand.

5. To obtain one estimate of the value of the integral of f (x) over the interval [a, b]

given n total random points in the rectangle containing the graph of f (x) and u

random points under the graph of f (x), calculate the percentage of points under

the curve, then multiply the area of the rectangle by the percentage: estimate =

× w × h.

To determine a conﬁdence interval for the approximation:

1. Create a list to store nvals results, where nvals is greater than 30.

2. In a loop, execute one of the functions described below and store the result in the

list.

For a uniform distribution, the probability of picking any value k from an interval [a,b] is the same for

all values k1, k2, ..., kn in the interval.

3. Execute the zInterval command by selecting the menu item Document Tools -

Statistics - Conﬁdence Interval - z interval ... Use the function stDevPop(list)

to specify the required population standard deviation.

2.1.1 Functions for ”Hit or Miss” Monte Carlo Integration

Based on the above description and observations, implementing ”hit or miss” Monte

Carlo integration with TI-Nspire functions is straight-forward. The attached TI-Nspire

document montecarlo.tns contains implementations of the following three functions

along with examples of using them to approximate integrals.

1. mcxypoints() - returns a matrix of (x, y) points with x and y values under the

graph of a function f (x) in rows 1 and 2, and x and y on or above the graph in

rows 3 and 4. The pseudocode for this function is presented in algorithm 1.

2. mcintegral() - returns the approximate value of the integral of a function f (x)

as the area under the graph of f (x). The pseudocode for this function is presented

in algorithm 2.

3. mcintxyz() - returns the approximate value of the integral of a function f (x, y, z)

as the hyper-volume under the graph of f (x, y, z). This function demonstrates the

ease of extending the Monte Carlo Integration technique for multi-variable func-

tions. The pseudocode for this function is presented in algorithm 3.

2.1.2 Examples of ”Hit or Miss” Monte Carlo Integration

The following examples show how to use the functions described above to ﬁnd the ap-

proximate value of an integral with Monte Carlo Integration.

Example 1.

a. Generate 1000 random points with the function mcxypoints() for approximat-

ing the value of

(x − 2)dx.

b. Estimate the value of the integral using the percentage of points under the graph.

c. Plot the function, the containing rectangle, and the points over and under the graph

of the function.

Solution:

First, deﬁne the function and ﬁnd the actual value of the integral in a calculator page

for comparison with the estimated value:

Next, determine the maximimum value of the function in the interval [2, 5]: since the

function is increasing on the interval, the maximum value occurs at the right-hand side

of the interval, and is f 2(5) = 5 − 2 = 3.

Finally, call the function mcxypoints() to return a matrix containing the random

points:

Extract the x, y values from the matrix as lists:

Note: The underline ( ) entries in the matrix and lists are TI-Nspire Void values, which

are ignored when the lists are graphed and by most TI-Nspire list functions.

Calculate the percentage of random points under the graph then multiply the area of

the rectangle by the percentage to estimate the area under the graph. Use the TI-Nspire

function count() to determine the number of points, then divide the count by 1000:

The lists of random points are plotted in a graphs page as scatter plots. Figure 3 shows

the graph of the containing rectangle, the function, and the random points.

Discussion

Figure 3: Random Points for Approximating the Integral of f 2(x) = x − 2

The value of the integral is 4.5 and the estimated value from Monte Carlo integra-

tion is 4.383. The error of the estimate is 4.5 − 4.383 = 0.117, which is a fairly large

error. However, the estimate is a single, random point estimate. Because each estimate

is random, each estimate will differ: some will be larger than the actual value of the

integral and some will be smaller. There are two ways to determine a more accurate

estimate: either increase the sample size, or take multiple samples and calculate the

mean (average) of areas from the multiple samples.

The ”hit or miss” Monte Carlo integration method is a probabilistic experiment where a

random sample of n trials is taken, with each trial consisting of random x, y coordinates

drawn from a uniform probability distribution. The result of the experiment is a sample

proportion ˆp: the percentage of the random points under the graph of a function. The

calculated area for each random sample is also a random variable.

By the strong law of large numbers, as the sample size n approaches inﬁnity, the

average proportion of the observations approaches the actual proportion of observa-

tions under the graph. Thus, increasing the number of trials for a sample theoretically

increases the accuracy of the calculated area under the graph. Example 2 demonstrates

using a large number of random points to approximate the value of an integral.

The second technique (calculating the mean of multiple samples) is justiﬁed by the

central limit theorem. By the central limit theorem, a distribution of multiple random

variables is normally distributed regardless of the underlying probability distribution

function from which the random samples were drawn, provided the number of samples

is large enough (greater than about 30 samples). The mean of the normally distributed

areas approaches the population mean (the mean area under the graph of the function).

As a result, multiple small samples can be used to approximate the value of an integral

by ﬁnding the mean of the distribution of sample areas instead of ﬁnding the area from

a single large sample. In addition, since the sample areas are normally distributed, a

conﬁdence interval for the range of values can also be calculated. Example 3 illustrates

this technique.

Example 2. Find the approximate value of

(x − 2)dx using ”hit or miss” Monte

Carlo integration with 10, 000 random points.

Solution: Run the function mcintegral() to ﬁnd the approximate value of the in-

tegral:

Discussion: The approximate value with a sample of 10, 000 random points is 4.5162

which is slightly higher than the actual value of 4.5. This is a much better approxima-

tion than the approximation with a sample of 1, 000 points. The error in the approxi-

mation is 0.0162.

Example 3. Find the approximate value of

(x − 2)dx using ”hit or miss” Monte

Carlo integration by ﬁnding the mean area of 100 samples, each consisting of 100 ran-

dom points.

Solution: Create a list to hold 100 areas, then run the function mcintegral() 100

times to ﬁnd the approximate value of the integral. Store the returned areas in the list,

then ﬁnd the mean of the list of areas:

Discussion: The approximate value of the integral from 100 samples of 100 random

points is 4.4811 which is slightly lower than the actual value of 4.5. This is a much

better approximation than the approximation with a sample of 1, 000 points. The error

in the approximation is 0.01889.

The returned list of approximate areas can be added to a Lists & Spreadsheet ap-

plication, a conﬁdence interval can be calculated, and the probability distribution of

the random values can be generated, as showm in Figure 4. The 95 percent conﬁdence

interval for the value of the area is [4.399, 4.563]. This conﬁdence interval means that

if the experiment is repeated 100 times, then about 95 of the times the actual value of

the mean area will be somewhere between 4.4 and 4.6.

Figure 4: Statistics for Approximating the Integral of f 2(x) = x − 2

Example 4. Use ”hit or miss” Monte Carlo Integration to ﬁnd the approximate value

of the integral

2.5

0.5



1+2cosh(2x)(ln(x))



dx.

Solution: The maximum value of the integrand in the interval of integration is needed

to ﬁnd the approximate value of the integral. The simplest way to determine this is by

graphing the integrand, selecting the graph, then picking the Analyze Graph - Maxi-

mum right-click popup menu item. First, deﬁne the function in a calculator page and

evaluate its integral to compare the actual value with the approximate value:

Add a graph page, graph the function, and ﬁnd its maximum value. Figure 5 shows

the graph of f 3(x) with the maximum value and the value of the integral from x = 0.5

to x = 2.5.

Now calculate two approximate values of the integral. Find the ﬁrst approximate

value as a single point estimate using 10,000 random points:

Figure 5: The Graph of f 3(x)

Find the second approximate value as the mean of 100 areas, each approximated using

100 points:

Calculate a 95 percent conﬁdence interval for the mean area:

Discussion: The error in the ﬁrst approximate value is 0.02. The error in the sec-

ond approximate value is 0.012 and the 95 percent conﬁdence interval is [0.704, 0.743].

Example 5.

Use ”hit or miss” Monte Carlo Integration to ﬁnd the approximate value

of the integral



−

(

)



dx dy dz.

Solution:

First, deﬁne the integrand in a calculator page and calculate the value of the integral:

Determine the maximum value of f xyz() on the interval [0, 1] × [0, 1] × [0, 1]: since

this is a function of three variables, its graph is in R

, so ﬁnding its maximum value

graphically is not feasible. Examining the function reveals it is simply an exponential

function of the form e

whose value when k = 0 equals 1.0. When k < 0, e

is a de-

creasing function that approaches 0 as k increases. Therefore, the maximum value of

the function over the interval of integration is 1.0.

Obtain a single point estimate of the value of the integral with 10, 000 random points

by running the three-variable implementation of the Monte Carlo integration function:

Find a second approximate value as the mean of 100 areas, each approximated us-

ing 100 points:

Calculate a 95 percent conﬁdence interval for the mean area:

Example from [Farlow], problem 3., page 345

Discussion: The error of the ﬁrst approximation is 0.0034 and is 0.0032 for the sec-

ond approximation. The actual value of the integral, 0.4165, is in the conﬁdence in-

terval, [0.4105, 0.4289]. Both techniques for approximating the value of the multi-

dimensional integral yield good results.

The implementation of the function mcintxyz() used for this example illustrates how

to easily extend the basic Monte Carlo Integration technique to approximate the value

of higher-dimension integrals. For a function of one variable, the percentage of a two-

dimensional area is calculated; for a function of two variables, the percentage of a

three-dimensional volume is calculated; for a function of three variables, the percent-

age of a four-dimensional hyper-volume is calculated; etc.

3 The Sample Mean Monte Carlo Integration Method

The sample mean Monte Carlo integration method is a simpler method than the ”hit or

miss” Monte Carlo integration method. It is also computationally more efﬁcient. The

method is derived from the deﬁnition of the expected value of a function of a random

variable.

The Law of the Unconscious Statistician asserts that the expected value of a contin-

uous function f (X), where X is a random variable from a probability density function

pd f (x), is given by

E[ f (X

)] = µ =

( f (x))(pd f (x)) dx

where R is the region of integration[wikipedia].

This deﬁnition of E[ f (X

)] is valid for multi-dimensional functions as well as for func-

tions of a single variable, as long as each random variable is from the same probability

density function.

The pdf for a random variable from a uniform probability density function over an

interval [a, b] is

b−a

. The formula for the expected value of f (X

) when X

is from a

uniform pdf is thus

E[Y

] = E[ f (X

)] = µ =



f (x))



b − a



Multiplying the expected value of the function of a random variable, E[ f (X

)] by (b− a)

results in a formula for ﬁnding an approximation to f (x) [Chihara]:

(b − a) · E( f (X

)) = µ = (b − a)



f (x))



b − a



dx =

( f (x)) dx

By the strong law of large numbers,

lim

N→∞

= µ =

( f (x)) dx

where

+ · · · +Y

) and Y

= (b − a) f (X

Thus, for large N,

≈

( f (x)) dx.

Based on the above derivation, the simple formula to compute

, the approximation

( f (x)) dx is

( f (x)) dx ≈

∑

(b − a) f (X

)

(b − a)

∑

f (X

)

where X

is a random variable from a uniform probability distribution over the interval

[a, b].

Although the derivation of the ﬁnal formula appears complicated, there is a simple

geometric interpretation of the formula: the value of each Y

is simply the area of a

rectangle of width b − a and height f (X

) and the ﬁnal result of the computation is the

mean (average) of the N rectangle areas.

Recall the First Mean Value Theorem for Integrals:

for some value c in [a, b],

( f (x)) dx = f (c)(b − a)

where f (c) is the average value of f (x) on [a, b] and f (c)(b − a) is the area of a rect-

angle of width b − a and height f (c). Geometrically, the sample mean Monte Carlo

integration method approximates f (c)(b − a) by averaging the areas of N randomly-

generated rectangles of width b − a and height f (X

3.1 Sample Mean Monte Carlo Integration with TI-Nspire

Since sample mean Monte Carlo integration does not require ﬁnding the maximum

value of the integrand on the interval of integration, the sample mean technique is

easier to perform than the ”hit or miss” technique:

1. Deﬁne the integrand in a calculator page.

2. Run one of the user-deﬁned functions described in the next section.

To determine a conﬁdence interval for the approximation:

1. Create a list to store nvals results, where nvals is greater than 30.

2. In a loop, execute one of the functions described below and store the result in the

list.

3. Execute the zInterval command by selecting the menu item Document Tools -

Statistics - Conﬁdence Interval - z interval ... Use the function stDevPop(list)

to specify the required population standard deviation.

3.1.1 Functions for Sample Mean Monte Carlo Integration

The attached TI-Nspire document montecarlo.tns contains implementations of the fol-

lowing two functions along with examples of using them to approximate integrals.

1. smmcintegral() - calculates and returns an approximation to a single-variable

integral using the computational formula

(b−a)

∑

f (X

). The pseudocode for

this function is presented in algorithm 4.

2. smmcintxyz() - calculates and returns an approximation to a function of three

variables by extending the above computational formula. The pseudocode for

this function is presented in algorithm 5.

3.1.2 Examples of Sample Mean Monte Carlo Integration

The following examples demonstrate using sample mean Monte Carlo integration to

approximate the same functions approximated above with ”hit or miss” Monte Carlo

integration. The results from ”hit or miss” Monte Carlo integration for each example

are presented for comparison with the results from sample mean Monte Carlo integra-

tion.

Example 6. Find the approximate value of

(x − 2)dx using sample mean Monte

Carlo integration with 10, 000 random points.

Solution: Run the function smmcintegral() to ﬁnd the approximate value of the

integral:

Discussion: The approximate value with a sample of 10, 000 random points is 4.4973

which is slightly less than the actual value of 4.5. The error in the approximation is

0.0027. From Example 2, the value with ”hit or miss” integration was calculated as

4.5162 with an error of 0.0162.

Example 7. Find the approximate value of

(x − 2)dx using sample mean Monte

Carlo integration by ﬁnding the mean area of 100 samples, each consisting of 100 ran-

dom points, then ﬁnd a 95 percent conﬁdence interval for the mean value.

Solution: Create a list to hold 100 areas, then run the function smmcintegral() 100

times to ﬁnd the approximate value of the integral. Store the returned areas in the list,

then ﬁnd the mean of the list of areas:

A conﬁdence interval is found by running the zInterval command:

Discussion: The approximate value of the integral from 100 samples of 100 random

points is 4.51462 which is greaater than the actual value of 4.5. The error in the approx-

imation is 0.01462 which is a larger error than the error from a single approximation

with 10000 random values. From Example 3, the value from ”hit or miss” integration

was calculated as 4.4811 with an error of 0.01899.

The 95 percent conﬁdence interval for the value of the area is [4.465, 4.5643]. This

conﬁdence interval means that if the experiment is repeated 100 times, then about 95 of

the times the actual value of the mean area will be somewhere between 4.47 and 4.56.

In Example 3, the conﬁdence interval for ”hit or miss” integration was [4.399, 4.563].

Example 8. Use sample mean Monte Carlo Integration to ﬁnd the approximate value

of the integral

2.5

0.5



1+2cosh(2x)(ln(x))



dx.

Solution: Run the function smmcintegral() to ﬁnd the approximate value of the

integral:

Discussion: The approximate value with a sample of 10, 000 random points is 0.73889

which is greater than the actual value of 0.735118. The error in the approximation is

about 0.0038. The ”hit or miss” value in Example 4 was 0.755 and the error was 0.02.

Example 9. Find the approximate value of

2.5

0.5



1+2cosh(2x)(ln(x))



dx using sample

mean Monte Carlo integration by ﬁnding the mean area of 100 samples, each consist-

ing of 100 random points, then ﬁnd a 95 percent conﬁdence interval for the mean value.

Solution: Create a list to hold 100 areas, then run the function smmcintegral() 100

times to ﬁnd the approximate value of the integral. Store the returned areas in the list,

then ﬁnd the mean of the list of areas:

A conﬁdence interval is found by running the zInterval command:

Discussion: The approximate value of the integral from 100 samples of 100 random

points is 0.73427 which is less than than the actual value of 73512. The error in the ap-

proximation is 0.00085 which is less than the error from a single approximation with

10000 random values. The error in the approximation is about 0.0038. The ”hit or

miss” value in Example 4 was 0.7234 and the error was 0.012.

The 95 percent conﬁdence interval for the value of the area is [0.72187,0.74668]. This

conﬁdence interval means that if the experiment is repeated 100 times, then about 95

of the times the actual value of the mean area will be somewhere between 0.722 and

0.747. The conﬁdence interval from ”hit or miss” integration was [0.7039, 0.7429].

Example 10.

Use sample mean Monte Carlo Integration to ﬁnd the approximate

value of the integral



−

(

)



dx dy dz.

Example from [Farlow], problem 3., page 345

Solution: Run the function smmcintxyz() to ﬁnd the approximate value of the in-

tegral:

Discussion: The approximate value with a sample of 10, 000 random points is 0.416889

which is greater than the actual value of 0.416538. The error in the approximation is

about 0.000351. The approximate value from the ”hit or miss” method in example 5

was 0.4199 with an error of .0034.

Example 11. Find the approximate value of



−

(

)



dx dy dz using

sample mean Monte Carlo integration by ﬁnding the mean area of 100 samples, each

consisting of 100 random points, then ﬁnd a 95 percent conﬁdence interval for the

mean value.

Solution: Create a list to hold 100 areas, then run the function smmcintxyz() 100

times to ﬁnd the approximate value of the integral. Store the returned areas in the list,

then ﬁnd the mean of the list of areas:

A conﬁdence interval is found by running the zInterval command:

Discussion: The approximate value calculated as the mean of 100 samples of 100

random points is 0.418333 which is greater than the actual value of 0.416538. The

error in the approximation is about 0.0018. From the ”hit or miss” method in example

5, the approximate value was 0.4197 and the error was .0032.

The 95 percent conﬁdence interval for the value of the integral is [0.4151,0.4216]. The

conﬁdence interval from the ”hit or miss” method in example 5 was [0.4197, 0.4289].

4 Observations about Monte Carlo Integration

4.1 Advantages and Disadvantages

Advantages:

• Can provide approximate values of integrals when other numerical techniques

can not be used.

• Easy to implement, extend, and use.

• Methods execute in linear time regardless of the number of iterations.

• Useful for approximating multi-variable integrals.

Disadvantages:

• Not as accurate as numerical methods, especially for single-variable integrals.

• Many iterations required to achieve accurate results.

See [Scratchapixel] for additional details about Monte Carlo methods and Monte Carlo

integration, including techniques for increasing the accuracy of Monte Carlo integra-

tion.

4.2 Comparing the Two Methods

Examining the errors from the two methods, the sample mean Monte Carlo method

yields better approximations than does the ”hit or miss” Monte Carlo method. In addi-

tion, since the sample mean method does not require determining the maximum value

of the integrand over the interval of integration, the sample mean method is easier to

use.

For a detailed comparison of the methods and an analysis of errors vs sample size

and processing time vs sample size for both methods, see the slides at the link https:

//sites.oxy.edu/ron/math/400/15/comps/Klenha.pdf.

5 Algorithms

5.1 Algorithms for ”Hit or Miss” Monte Carlo Integration

Algorithm 1 mcxypoints - Generate Matrix of Random Point Values

Input:

f is the name of f(x), the integrand

f max is the max value of f(x) on [a,b]

a is the lower bound of the integral

b is the upper bound of the integral

npts is the number of points to generate

Output: xy matrix of random points

function MCXYPOINTS( f , f max, a, b, npts)

xy ← newMat(4, npts)

interval ← b − a

for j = 1, j ≤ npts do

x ← a + rand()· interval  random x in [a,b]

y ← rand() · f max  random y in [0,fmax]

f o f x ← # f (x)  f (randomx)

if y < f o f x then  x, y is below the graph

xy[1, j] ← x

xy[2, j] ← y

xy[3, j] ← Void

xy[4, j] ← Void

else  x, y is above or on the graph

xy[1, j] ← Void

xy[2, j] ← Void

xy[3, j] ← x

xy[4, j] ← y

end if

end for

return xy

end function

Algorithm 2 mcintegral - Calculate the approximate value of an integral.

Input:

f is the name of f(x), the integrand

f max the max value of f(x) on [a,b]

a is the lower bound of the integral

b is the upper bound of the integral

npts is the number of points to generate

Output: approximate value of the integral of f(x) over [a,b]

function MCINTEGRAL( f , f max, a, b, npts)

under ← 0

interval ← b − a

for j = 1, j ≤ npts do

x ← a + rand()· interval  random x in [a,b]

y ← rand() · f max  random y in [0,fmax]

f o f x ← # f (x)  f (randomx)

if y < f o f x then  x, y is below the graph

under ← under + 1

end if

end for

return

under

npts

· f max · interval

end function

Algorithm 3 mcintxyz - Calculate the approximate value of an integral of 3 variables.

Input:

f is the name of f(x), the integrand

f max the max value of f(x) on [a,b]

xa, xb is the lower and upper x bounds of the integral

ya, yb is the lower and upper y bounds of the integral

za, zb is the lower and upper z bounds of the integral

npts is the number of points to generate

Output: approximate value of the integral of f(x,y,z) over [xa, xb] × [ya, yb] × [za, zb]

function MCINTXYZ( f , f max, xa, xb,ya, yb, za, zb, npts)

under ← 0

xinterval ← xb − xa

yinterval ← yb − ya

zinterval ← zb − za

for j = 1, j ≤ npts do

x ← xa + rand()· xinterval  random x in [xa,xb]

y ← ya + r and() · yinterval  random y in [ya,yb]

z ← za + rand() · zinterval  random y in [za,zb]

w ← f max · rand()  random w

f val ← # f (x, y, z)  value of f(x,y,z)

if w < f val then  f (x, y, z) is below the graph

under ← under + 1

end if

end for

return

under

npts

· f max · xinterval · yinterval · zinterval

end function

5.2 Algorithms for Sample Mean Monte Carlo Integration

Algorithm 4 smmcintegral - Calculate the approximate value of an integral.

Input:

f is the name of f(x), the integrand

a is the lower bound of the integral

b is the upper bound of the integral

npts is the number of points to generate

Output: approximate value of the integral of f(x) over [a,b]

function SMMCINTEGRAL( f , a, b, npts)

sumvals ← 0

interval ← b − a

for j = 1, j ≤ npts do

x ← a + rand()· interval  random x in [a,b]

sumvals ← sumvals + # f (x)  add f (randomx) to sum

end for

return

interval·sumvals

npts

end function

Algorithm 5 smmcintxyz - Calculate the approximate value of an integral of 3 vari-

ables.

Input:

f is the name of f(x), the integrand

xa, xb is the lower and upper x bounds of the integral

ya, yb is the lower and upper y bounds of the integral

za, zb is the lower and upper z bounds of the integral

npts is the number of points to generate

Output: approximate value of the integral of f(x,y,z) over [xa, xb] × [ya, yb] × [za, zb]

function SMMCINTXYZ( f , xa, xb, ya, yb, za, zb, npts)

sumvals ← 0

xinterval ← xb − xa

yinterval ← yb − ya

zinterval ← zb − za

interval ← xinterval · yinterval · zinterval

for j = 1, j ≤ npts do

x ← xa + rand()· xinterval  random x in [xa,xb]

y ← ya + r and() · yinterval  random y in [ya,yb]

z ← za + r and() · zinterval  random y in [za,zb]

sumvals ← sumvals + # f (x, y, z)  add f(x,y,z) to sum

end for

return

sumvals·interval

npts

end function

References

[Chihara] Chihara, Laura, and Hesterberg, Tim, Mathematical Statistics with Resam-

pling and R, John Wiley & Sons, Inc., Hoboken, NJ, 2011, pp. 341-345

[Farlow] Farlow, Stanley J., Partial Differential Equations for Scientists and Engi-

neers, Dover Publications, Inc., New York, NY, 1993, pp. 340-345

[Napolitano] Napolitano, Jim, A Mathematica Primer for Physicists, CRC Press, Tay-

lor & Francis Group, Boca Raton, FL, 2018, pp. 136-138

[Scratchapixel] Monte Carlo Methods in Practice. Scratchapixel 2.0,

www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/

monte-carlo-methods-in-practice/monte-carlo-integration.

Accessed 14 Oct 2020

[wikipedia] Expected Value. wikipedia.org,

https://en.wikipedia.org/wiki/Expected_value.

Accessed 29 Oct 2020