Lecture : 6

Computer Simulations and Monte Carlo

Methods
 Problem
definition

PLANNING
Strategic/tactical
planning Model Scoping

Model building
Data acquisition
Model Modeling

Validation
VERIFICATION

Experiments

Analysis
APPLICATION
Fig.1 Overview of the simulation methodology
Fundamentals of Performance evaluation of computer and
Telecommunication Systems-Obaidat & Boudriga, John wiley, 2010 Implementation/
documentation
 Actual World Study of Simulation
System
Model
under study

Simulation
experiment

Analysis of
simulation

Result

Changed system

Fig.2 Overall simulation process, ibid

 Formulation
of Problem Design

Objective and Final Analysis and

project plan setting production runs
Yes
Mode Collection More
structure of data runs

no
Code Study of document
no and reporting
no results

If verified

If validated Implementation
no Yes

Fig.3 The simulation modeling process. (ibid)

 Ways to Study System
System

Experiment w/ Experiment w/
actual system model of system

Physical Mathematical
Model Model

Analytical Simulation
Model
Model

Fig.4
Focus of class
5
Probability and Statistics for Computer Scientists-Michael Baron, CRC,2007.

Computer simulations refer to a generation of a process by writing a suitable

computer program and observing its results. Monte Carlo methods are those
based on computer simulations involving random numbers.

The main purpose of simulations is computing such quantities of interest whose

direct computation is complicated, risky, consuming , expensive, or impossible. For
example, suppose a complex device or machine is to be built and launched.
Before it happens, its performance is simulated allowing experts to evaluate its
adequacy and associated risks carefully and safely. For example, one surely prefers
to evaluate safety and reliability of a new module of a space station by means of
computer simulations rather than during the actual mission.

Monte Carlo methods are mostly used for the computation of probabilities,
expected values, and other distribution characteristics. Recall that probability can
be defined as long-run proportion. With the help of random number generators,
computer can actually simulate long run. Then, probability can be estimated by a
mere computation of the associated frequency. Similarly, one can estimate
expectations, variances, and other distribution characteristics from a long run of
simulated random variables.
 What is Monte Carlo ?
Monte Carlo is an administrative area of the principality of
Monaco.

Famous for its casinos !

Monte Carlo is a (large) class of numerical methods used to solve

integrals and differential equations using sampling and
probabilistic criteria.
7
 Introduction
• Monte Carlo methods are stochastic
techniques.
• Monte Carlo method is very general.
• We can find MC methods used in
everything from economics to nuclear
physics to regulating the flow of traffic.
 Introduction
• A Monte Carlo method can be loosely
described as a statistical method used in
simulation of data.
• And a simulation is defined to be a method
that utilizes sequences of random numbers as
data.
 Introduction (cont.)

• The Monte Carlo method provides approximate

solutions to a variety of mathematical
problems by performing statistical sampling
experiments on a computer.
• The method applies to problems with no
probabilistic content as well as to those with
inherent probabilistic structure.
 Monte Carlo Simulations
A Monte Carlo method is a computational
algorithm that relies on repeated random
sampling to compute its results.
Nicholas Constantine Metropolis

Worked at Los Alamos on the Manhatten

project .Worked with Enrico Fermi
(developed quantum physics) and Edward
Teller (father of the hydrogen bomb) on the
first nuclear reactors.

John von Neumann

The term Monte Carlo was coined in the 1940s by physicists working on
nuclear weapon projects in the Los Alamos National Laboratory.
Andrei A. Markov 1856-1922
Russian
mathematicians

A. N. Kolmogorov (12.4.1903-20.10.1987)
• Monte Carlo methods are often used when
simulating physical and mathematical systems.

•Because of their reliance on repeated

computation and random or pseudo-random
numbers, Monte Carlo methods are most
suited to calculation by a computer. Monte
Carlo methods tend to be used when it is
infeasible or impossible to compute an exact
result with a deterministic algorithm.
 The simplest Monte Carlo method
Finding the value of π (“shooting darts”)
π/4 is equal to the area of a
circle of diameter 1.

1 4
  dx dy
x 2  y 2 1

Integral solved
 with Monte Carlo
1

Details of the Method

Randomly select a large number of points inside the square
N _ inside_ circle (Exact when N_total  ∞)
  4*
N _ total
15

Monte Carlo Example:
Estimating 
If you are a very poor dart player, it is easy to imagine throwing
darts randomly at the above figure, and it should be apparent that of
the total number of darts that hit within the square, the number of
darts that hit the shaded part (circle quadrant) is proportional to the
area of that part. In other words,
If you remember your geometry, it's easy to show that
 (x, y)

x = (random#)
y = (random#)
distance = sqrt (x^2 + y^2)
if distance.from.origin (less.than.or.equal.to) 1.0
let hits = hits + 1.0
 Cont....
area.of .circle #.of .dots.inside.circle

area.of .square total.number.of .dots

Hit and miss algorithm

Generate two sequences of N of PRN :: Ri,,Rj
Xi=-1+2R i
Yj=-1+2R j
Start from s=zero
If (X²+Y²<1) s=s+1
# of dots inside circle=s
total number of dots=N
  4*S / N
 Monte Carlo Integration
♥ Hit and miss method
♥ Sample mean method
 Hit and Miss method
a, b  R
b

I   f ( x) dx
a

♦Generate two sequence of N of PRN (Ri,R j) i& j=1,2,….,N

♦0 ≤f(x) ≤Ymax ,for X Є (a,b)

♦ Xi=a+Ri (b-a)
♦ Yi=Ymax R j
♦ start from s=0
♦ if Yj<f(x) s=s+1
♦ I=Ymax(b-a) S/N
program Hit_miss
implicit none
real(8),dimension(:),allocatable::R
real(8)::X,Y,m,a,b,integ
integer::i,N,s
read*,a !the lower value of the interval
read*,b !the upper value of the interval
M=f(b)
if (dabs(f(a))>dabs(f(b))) M=f(a)
read*,N !the number of random numbers
allocate(R(n))
s=0
do i=1,N
call random_number(R(n))
call random_number(R(n+1))
X=a+R(n)*(b-a)
Y=M*R(n+1)
if (y<f(x)) s=s+1
end do
INTEG=M*(b-a)*s/N
print*,integ
contains
real(8) function F(X)
real(8),intent(in)::X
F=2.0*X+1.0
end function F
end program Hit_miss
 Cont….
If there is a part of the function under X-axis
b

s area.under.curve  f ( x)dx
  a

N Total.area ( M  R)(b  a)
b
s
 f ( x)dx  ( M  R)(b  a)
a N
This must now be corrected for the fact that the line y = -R has
been used as the baseline for the integral instead of the line
y = 0. This is accomplished by subtracting the rectangular
area R(b-a). b
s
then  f ( x)dx  ( M  R)(b  a)  R(b  a)
a N
 Types of distribution

Uniform distribution Gaussian or normal distribution

 MONTE CARLO INTEGRATION
AND
VARIATION REDUCTION METHOD
(for random quadrature)
𝒃
𝑰 = න 𝒈 𝒙 𝒅𝒙
𝒂

Sample Mean Method

We may write 𝑰 as:

𝒃 𝒃
𝒈 𝒙
𝑰 = න 𝒈 𝒙 𝒅𝒙 = න 𝒇 𝒙 𝒅𝒙
𝒂 𝒂 𝒇 𝒙
for some pdf 𝒇(𝒙) with support [𝒂, 𝒃].
Now if 𝑿𝟏 , 𝑿𝟐 , 𝑿𝟑 , … … 𝑿𝒏 are iid uniformly distributed on the interval
𝒂, 𝒃 ,then
𝒃
𝑱 = න 𝒉 𝒙 𝒅𝒙 ≈ 𝒃 − 𝒂 < 𝒉 >
𝒂
where
𝒏
𝟏
< 𝒉 > = ෍ 𝒉 𝒙𝒊 .
𝒏
𝒊=𝟏
Based on the above, if we choose
𝟏
𝒇 𝒙 =
𝒃−𝒂
𝒈(𝒙)
Then = 𝒃 − 𝒂 𝒈(𝒙) and 𝑰 could be written as:
𝒇(𝒙)
(𝒃−𝒂) 𝒏
𝑰෠ = σ𝒊=𝟏 𝒈(𝒙𝒊 )
𝒏
 (𝑏−𝑎) 𝑛
𝐼መ = σ𝑖=1 𝑔 𝑥𝑖
𝑛
Using the foregoing , 𝐼መ approximates 𝐼.
መ
The error in the approximation can be calculated from the variance of 𝐼.

𝟏
ഥ = (𝑿𝟏 + 𝑿𝟐 + ⋯ … . +𝑿𝒏 ) and the distribution of 𝑿𝟏 , 𝑿𝟐 , … . , 𝑿𝒏 has a mean 𝝁 and standard deviation 𝝈, we
Since 𝑿 𝒏
observe that
𝟏 𝟏
𝑬 𝑿ഥ = 𝑬𝑿𝟏 + 𝑬𝑿𝟐 + ⋯ . +𝑬𝑿𝒏 = 𝒏𝝁 = 𝝁
𝒏 𝒏
And
𝟏 𝟏 𝝈𝟐
𝒗𝒂𝒓 𝑿 ഥ = 𝒗𝒂𝒓𝑿 𝟏 + 𝒗𝒂𝒓𝑿 𝟐 + ⋯ . +𝒗𝒂𝒓𝑿 𝒏 = 𝒏𝝈 𝟐
=
𝒏𝟐 𝒏𝟐 𝒏

In view of the above,

𝑛 F(x)
(𝑏 − 𝑎)
𝑣𝑎𝑟 𝐼መ = 𝑣𝑎𝑟
𝑛
෍ 𝑔 𝑥𝑖 = (𝑏 − 𝑎)2 𝑣𝑎𝑟 𝑔ҧ f
𝑖=1
(𝑏 − 𝑎)2
= 𝑣𝑎𝑟 𝑔
𝑛
1
where 𝑣𝑎𝑟 𝑔 =< 𝑔2 > −< 𝑔 >2 ≈ 𝑛 σ𝑛𝑖=1 𝑔2 𝑥𝑖 −
1 𝑛 2
𝑛
σ𝑖=1 𝑔(𝑥𝑖 ) x
a b
 Sample Mean MC algorithm
♠ Generate sequence of N of PRN : Ri
♠Compute Xi=a+Ri (b-a)
♠ compute f(Xi) , i=1,2,3,….,N
♠ use ˆ (b  a ) N
I   f (x ) i
N i 1

♠♠♠ note:: if f(x) is not square integrable ,then the MC

Estimate Î will still converge to the true value, but
The error estimate becomes unreliable.
program Sampled_Mean
implicit none
real(8),dimension(:),allocatable::R
real(8)::a,b,sum,integ,X
integer::i,N
read*,a ! Lower value
read*,b ! upper value
read*,n ! number of random points
allocate(R(n))
sum=0.0d0
do i=1,n
call random_number(R(n))
call random_number(R(n+1))
X=a+R(n)*(b-a)
sum=sum+f(x)
int=((b-a)/N)*sum
end do
write(*,*) "integ=",integ
contains
real(8) function F(X)
real(8),intent(in)::X
F=2*X+1.0d0
end function F
end program Sampled_Mean
 Numerical Example of Monte Carlo Integration
b
•Suppose interested in  sin( x )dx
0
– Simple problem with known solution
•Considerable variability in quality of solution for varying b
– Accuracy of numerical integration sensitive to integrand
and domain of integration

Integral estimates for varying n

n = 20 n = 200 n = 2000

b=
2.296 2.069 2.000
(ans.=2)
b = 2
0.847 0.091 0.0054
(ans.=0)
31
 Method of Generation of Random Numbers
Analysis of Computer and Communication Networks, Fayez Gebali, Springer,2008

Transforming Random variables

Mathematical packages usually have functions for generating random numbers
following the uniform and normal distributions ony. However, when we are
simulating communication systems, we need to generate network traffic that
follows other types of distributions. How can we do that? Well, we can do that
through transforming random variables, which is the subject of this section.

Continuous Case
Suppose we have a random variable 𝑿 with known pdf and cdf and we have
another random variable 𝒀 that is a function of 𝑿:
𝒀 = 𝒈(𝑿)
Here 𝑿 is named the source random variable and 𝒀 is named the target random
variable. We are interested in finding the pdf and cdf of 𝒀 when the pdf and cdf of
𝑿 are known. The probability that 𝑿 lies in the range 𝒙 and 𝒙 + 𝒅𝒙 is given by
𝑷 𝒙 ≤ 𝑿 ≤ 𝒙 + 𝒅𝒙 = 𝒇𝑿 𝒙 𝒅𝒙
But this probability must equal that 𝒀 lies in the range 𝒚 and 𝒚 + 𝒅𝒚. Thus we can
write
𝒇𝒀 𝒚 𝒅𝒚 = 𝒇𝑿 𝒙 𝒅𝒙
where 𝒇𝒀 (𝒚) is the pdf for the random variable 𝒀 and it was assumed that the
function 𝒈(𝒙) was monotonically increasing with 𝒙.
If 𝑔 was monotonically decreasing with 𝑥 , then we would have
𝑓𝑌 𝑦 𝑑𝑦 = −𝑓𝑋 𝑥 𝑑𝑥

The last two equations define the fundamental law of probabilities , which is given
by
𝑑𝑥
𝑓𝑌 𝑦 𝑑𝑦 = 𝑓𝑋 𝑥 𝑑𝑥 ⟹ 𝑓𝑌 𝑦 = 𝑓𝑋 𝑥
𝑑𝑦
since 𝑓𝑌 𝑦 and 𝑓𝑋 𝑥 are always positive.

In the discrete case the fundamental law of probability gives

𝑝𝑌 (𝑦) = 𝑝𝑋 (𝑥)
where 𝑝𝑋,𝑌 𝑥, 𝑦 is the pmf of the source and the target random variable
respectively.

# Given the random variable 𝑋 whose pdf has the form

2
𝑓𝑋 𝑥 = 𝑒 −𝑥 ; 𝑥 ≥ 0
Find the pdf for the random variable 𝑌 that is related to 𝑋 by the relation
𝑌 = 𝑋2
Sol.: We have
𝑑𝑥 1
𝑥=± 𝑦⟹ =±
𝑑𝑦 2 𝑦
Therefore
1 −𝑥 2
1
𝑓𝑌 𝑦 = 𝑒 ⟹ 𝑓𝑌 𝑦 = 𝑒 −𝑦 ; 𝑦 ≥ 0
2 𝑦 2 𝑦

# Suppose our source random variable is uniformly distributed and our target
random variable 𝑌 is to have a pdf given by
𝑓𝑌 𝑦 = 𝜆𝑒 −𝜆𝑦 ; 𝜆 > 0
Derive 𝑌 as a function of 𝑋.

Sol.: this example is fundamentally important since it shows how we can find the
transformation that allows us to obtain random numbers following a desired pdf
given the random numbers for the uniform distribution.

Using the relation

𝑑𝑥
𝑓𝑌 𝑦 = 𝑓𝑋 𝑥
𝑑𝑦
We get
𝑑𝑥
𝜆𝑒 −𝜆𝑦 = 𝑓𝑋 𝑥
𝑑𝑦
Assume the source random variable 𝑋 is confined to the range 0-1. The distribution

of 𝑓𝑋 𝑥 , shown in the Fig.1, suggests that 𝑓𝑋 𝑥

𝑓𝑋 𝑥 =1 and therefore we can write 1

𝑑𝑥
𝜆𝑒 −𝜆𝑦 =
𝑑𝑦
which on integration gives 0
𝑒 −𝜆𝑦 = 𝑥 0 1 𝑥
And we obtain the dependence of 𝑌 on 𝑋 as:
ln 𝑋
𝑌=𝑔 𝑋 =− ; 0 ≤ 𝑥 ≤ 1.
𝜆

Discrete Case
In the discrete case the fundamental law of probability gives
𝑝𝑌 𝑦 = 𝑝𝑋 𝑥
where 𝑝𝑋,𝑌 𝑥, 𝑦 corresponds to the pmf of the source / target random variable.
The values of 𝑌 are related to the values of 𝑋 by the equation
𝑌 = 𝑔(𝑋)
The procedure for deriving the pmf of 𝑌 given the pmf for 𝑋 is summarized as
follows:
1. For each value of 𝑥 obtain the corresponding value 𝑝 𝑋 = 𝑥 through
observation or modeling.
2. Calculate 𝑦 = 𝑔 𝑥 .
3. Associate 𝑦 with probability 𝑝 𝑋 = 𝑥 .

Generating Random Numbers

Random Computer Calculations???
• Compare this to the definition of an algorithm
(dictionary.com):
– algorithm
• n : a precise rule (or set of rules) specifying
how to solve some problem.
 Random Number
• What is random number ? Is 3 ?
– There is no such thing as single random number
• Random number
– A set of numbers that have nothing to do with the other
numbers in the sequence
• In a uniform distribution of random numbers in the
range [0,1] , every number has the same chance of
turning up.
– 0.00001 is just as likely as 0.5000
 Random v. Pseudo-random
• Random numbers have no defined sequence or
formulation. Thus, for any n random numbers, each
appears with equal probability.
• If we restrict ourselves to the set of 32-bit integers,
then our numbers will start to repeat after some very
large n. The numbers thus clump within this range
and around these integers.
• Due to this limitation, computer algorithms are
restricted to generating what we call pseudo-random
numbers.

June 19, 2020 OSU/CIS 541 39

 Monte-Carlo Methods
• 1953, Nicolaus Metropolis
• Monte Carlo method refers to any method
that makes use of random numbers
– Simulation of natural phenomena
– Simulation of experimental apparatus
– Numerical analysis
How to generate random numbers ?
• Use some chaotic system (Balls in a barrel – Lotto)
• Use a process that is inherently random
– Radioactive decay
– Thermal noise
– Cosmic ray arrival
• Tables of a few million random numbers
• Hooking up a random machine to a computer.
Pseudo Random number generators

• The closest random number generator that can be

obtained by computer algorithm.
• Usually a uniform distribution in the range [0,1]
• Most pseudo random number generators have two
things in common
– The use of large prime numbers
– The use of modulo arithmetic
• Algorithm generates integers between 0 and M, map
to zero and one.
X n  In / M
 An early example (John Von
Neumann,1946)
• To generate 10 digits of integer
– Start with one of 10 digits integers
– Square it and take middle 10 digits from answer
– Example: 57721566492 = 33317792380594909291
• The sequence appears to be random, but each number is determined
from the previous  not random.
• Serious problem : Small numbers (0 or 1) are lumped together, it can get
itself to a short loop. For example:
• 61002 = 37210000
• 21002 = 04410000
• 41002 = 16810000
• 51002 = 65610000
Anyone who considers arithmetical methods of producing random digits is, of
course, in a state of sin. (John Von Neumann,1951)
 Linear Congruential Method

• Lehmer, 1948
• Most typical so-called random number generator
• Algorithm : I n1  (aI n  c) mod(m)
– a,c >=0 , m > I0,a,c
• Advantage :
– Very fast
• Problem :
– Poor choice of the constants can lead to very poor sequence
– The relationship will repeat with a period no greater than m
(around m/4)
• C complier RAND_MAX : m = 32767
 RANDU Generator

• 1960’s IBM
• Algorithm

I n 1  (65539  I n ) mod(2 ) 31

• This generator was later found to have a

serious problem
1D and 2D Distribution of RANDU
Uniform Distribution
To generate a sequence of integer random numbers obeying the uniform
distribution using C programming , one invokes the function srand(seed). This
function returns an integer in the range 0 to RAND_MAX. When a continuous
random number is desired , drand is used. The function rand in MATLAB is
extensively used, in engineering & physics , to generate random numbers having
uniform distribution.

Inversion Method
To generate a number obeying a given distribution , we rely on the fundamental
law of probabilities discussed earlier. When 𝑋 has a uniform distribution over
the range 0 ≤ 𝑥 ≤ 1, we may use the following procedure for obtaining 𝑦.
We have
𝑑𝑥
= 𝑓(𝑦)
𝑑𝑦
Here 𝑓 𝑥 = 1 and 𝑓(𝑦) is the function describing the pdf of the target
distribution. Therefore integrating the foregoing
𝑦
‫׬‬0 𝑓 𝑧 𝑑𝑧 = 𝑥 ⟹ 𝐹 𝑦 = 𝑥 ⟹ 𝐹 −1 𝑥 = 𝑦.
The procedure then to obtain the random numbers corresponding to 𝑦 is to
generate 𝑥 according to any standard random number generator producing a
uniform distribution . Then use the above equation to provide us with 𝑦. Thus
the target random number value 𝑦 is obtained according to the following steps:
1. Obtain the source random number 𝑥 having a uniform distribution.
2. Lookup the value of 𝐹(𝑦) that is numerically equal to 𝑥 according to 𝐹 𝑦 = 𝑥.
3. Find the corresponding value of 𝑦 according to 𝐹 −1 𝑥 = y. If the inverse
function is not available , then lookup table is used and 𝑦 is chosen according to
the criterion 𝐹(𝑦) ≤ 𝑥 ≤ 𝐹(𝑦 + 𝜖), where 𝑦 + 𝜖 denotes the next entry in the
table.
Graphically the procedure is summarized in Fig.2
𝐹(𝑦)

Input 𝑥

𝑦
output 𝑦

Fig.2
 Sampling from a distribution
• We have a uniform distribution, how do we
sample from
– An exponential distribution?
– Other continuous distributions?
– A Poisson distribution?
 Exponential distribution
• Time between
radioactive decays.
• Distance a photon
travels before an
interaction. Probability density function (pdf)
• Mean 1/, variance
1/

(Wikipedia)

Cumulative distribution function (cdf)

 Sampling from the exponential
distribution
• Sample u from [0,1).
• Locate u on the y-axis
of the exponential
u
distribution cdf.
 Sampling from the exponential
distribution
• Sample u from [0,1).
• Locate u on y-axis.
• Find the x value that u
corresponds to u.

x
 Sampling from the exponential
distribution
• Sample u from [0,1).
• Locate u on y-axis.
• Find x = F-1(u).
u

x
 Acceptance-rejection

• To sample from a pdf

f(x):
1. Choose pdf g(x) and
constant c such that
c*g(x) ≥ f(x) for all x

www.mathworks.com
 Acceptance-rejection

• To sample from f(x):

1. Choose g(x) and c:
c*g(x) ≥ f(x) for all x.
2. Generate a random
number v from g(x).
3. Generate a random
number u from the
uniform distribution
on (0,1). v
 Acceptance-rejection

• To sample from f(x):

1. Choose g(x) and c:
c*g(x) ≥ f(x) for all x.
2. Generate v from g(x).
3. Generate u from
uniform distribution.
4. If c*u ≤ f(v)/g(v) accept
v
5. else reject v and go to
2.
 Acceptance-rejection

• To sample from f(x):

1. Choose g(x) and c:
c*g(x) ≥ f(x) for all x.
2. Generate v from g(x).
3. Generate uniform u.
4. If c*u ≤ f(v)/g(v) accept
v
5. else reject, go to 2.
• The average number
of iterations is c.