Monte Carlo Method (1)
What is a Monte Carlo simulation?
Monte Carlo simulation is a numerical method based on the extensive
use of random numbers for solving different problems.
Remark:
Monte Carlo simulation can be used to study not only stochastic
problems (e.g.,diffusion) but also non stochastic ones (Monte Carlo
integration).
A short history
Buffon’s problem of needle throwing:
A needle of length L is thrown at random onto a plane with
straight parallel lines which are separated by a distance d (d > L),
What is the probability P that the needle will intersect one of these
lines?
Georges Louis Leclerc comte de Buffon (1707 - 1788),Eminent
French naturalist,
Solution:
P = 2L/(pd)
Reference:
G,Comte de Buffon,Essai d ’arithmétique morale,Supplément à
l ’Histoire Naturelle,Vol,4,1777.
Remark:
Laplace suggested that the method described in Buffon’s problem
can be used as a stochastic method to calculate the value of p.
qh x
d
Condition for intersection:
h=Lsin(q)/2 > x
Geometric implication:
P = SI/(d/2?p/2)
SI = L/2
x
d/2
L/2
qp/2
x=Lsin(q)/2
I
Anecdote:
In 1901,Lazzerini (Italian mathematician) performed a simulation
by spinning round and dropping a needle 3407 times,He estimated
p to be 3.1415929 (accurate to the seventh number after the point!).
Exercise
Write a program to calculate p by using the method of Buffon with
d=1 and L=3/4.
Monte Carlo integration
x
d/2
L/2
qp/2
x=Lsin(q)/2
I
2/
0
s in2
p
qqdLS I
SI can be calculated from Pdp/4 (P is
the probability of finding a random
point in the area under the red curve).
Random number generation
General remark:
All the random numbers generated on a computer are in fact pseudo-
random numbers which are produced through a deterministic
algorithm!
Basic generator:
It produces uniformly distributed random numbers on [0,1],
Desirable properties of a random number generator:
It should generate sequences of random numbers which are
uniform,
uncorrelated,
with an extremely long period.
A widely used random-number generation method
Congruential method:
Xn = (aXn-1+c)mod(m) (1)
Properties of this sequence:
Xn < m
When a,c,X0 and m are integers,Eq.(1) generates integers between 0
and m-1,
Xn/m gives pseudo-random numbers distributed on (0,1).
Illustration:
With a=3,c=1,m=16 and a seed X0=2,we have the following
sequences,
2,7,6,3,10,15,14,11,2,7,…
The period of this sequence is 8.
Condition for a congruential generator to have a full period:
It can be shown that the congruential generator given in Eq.(1) has
a full period,i.e,m,if and only if
1) c and m have no common divisor (i.e.,no common divisor other
than 1);
2) a?1 mod(g) for every prime factor of m;
3) a?1 mod(4) if m is multiplier of 4.
It is obvious that m should be as large as possible! m=231-1 is often
used on a computer with 32 bits.
Statistical test of random number generators
Uniformity test,Break up the interval between 0 and 1 into a large
number of small bins and after generate a large number of random
numbers check for uniformity in the numbers of entries in each
bins (uniform histogram).
Parking lot test,Plot points in an m-dimensional space where the
m-coordinates of each point are determined by m-successive calls
to the random number generator,Then,look for regular structure
(resulting from some correlation).
Illustration:
good generator (parking lot test) bad generator(parking lot test)
1
100
1
100
How to generate distributions other than the
uniform one on (0,1)?
Transformation of random variables:
1) Uniform distribution on (a,b):
Let x be an uniform distribution on (0,1),the uniform distribution
on (a,b) is obtained from
y=a+(b-a)x
2) Random integers from 0 to N:
First divide the interval [0,1] into N+1 segments (each has a length
equal to D=1/(N+1)),M (0<M<N) is generated when
MD < x < (M+1)D
Code:
integer mran
mran=(N+1)*ran()
3) Random variables with an arbitrary distribution function
Let F(x) be the cumulative distribution function.
Remark:
F (x) is itself a random variable,Moreover,it can be shown that F(x)
is uniform random number on (0,1).
Algorithm:
1) generate a uniform random number,?1;
2)?2 = F-1(?1);
1 is a random number with the cumulative distribution F(x).
Application:
This method is most useful when F-1(x) can be obtained analytically.
Geometric illustration:
Example:
Histogram:
x
N(x)
0 1F(x) = x2
N(x)
0 1 x
N(x)
0 1
rescale the
abscissa
replot with
regular bins
1
2
x
f(x)
0
2x
xx d s sxF 2
0
2)(
Box-Müller method (generator of standard normal distribution):
Algorithm:
1) generate two uniform random numbers,?1,and?2 on (0,1);
2) n1 = (-2ln?1)1/2cos(2p?2)
n2 = (-2ln?1)1/2sin(2p?2)
give two independent random variables with standard Gaussian
distribution,i.e.,m=0 and s=1.
3) random variables with arbitrary m and s can be obtained from
n’1 = (n1 - m)/s
n’2 = (n2 - m)/s
Rejection techniques:
Basic rejection method:
Purpose:
Generate a sequence of random numbers which probability distribution
function is f(x),f(x) is bound by c and defined on a finite interval [a,b].
Algorithm:
1) generate one uniform random
number,?1,on (a,b) and one uniform
random number,?2,on (0,c);
2) accept?1 if?2? f(?1)/c;
3) repeat the above process.
Geometric meaning:
1 and?2 are the coordinates of random points over the rectangle,
When we keep only?1 for the points below f(x),?1 has the probability
distribution given by f(x).
c
a b
f(x)
x
Efficiency of rejection method:
The efficiency of a rejection is low when the surface area below f(x)
is small.
Illustration:
f1(x) can be sampled more efficiently by a rejection method that
f2(x).
c
a b
f2(x)
x
c
a b
f1(x)
x
Limitation of the method:
One must know the maximum of the distribution function and f(x)
should be defined on a finite interval.
What is a Monte Carlo simulation?
Monte Carlo simulation is a numerical method based on the extensive
use of random numbers for solving different problems.
Remark:
Monte Carlo simulation can be used to study not only stochastic
problems (e.g.,diffusion) but also non stochastic ones (Monte Carlo
integration).
A short history
Buffon’s problem of needle throwing:
A needle of length L is thrown at random onto a plane with
straight parallel lines which are separated by a distance d (d > L),
What is the probability P that the needle will intersect one of these
lines?
Georges Louis Leclerc comte de Buffon (1707 - 1788),Eminent
French naturalist,
Solution:
P = 2L/(pd)
Reference:
G,Comte de Buffon,Essai d ’arithmétique morale,Supplément à
l ’Histoire Naturelle,Vol,4,1777.
Remark:
Laplace suggested that the method described in Buffon’s problem
can be used as a stochastic method to calculate the value of p.
qh x
d
Condition for intersection:
h=Lsin(q)/2 > x
Geometric implication:
P = SI/(d/2?p/2)
SI = L/2
x
d/2
L/2
qp/2
x=Lsin(q)/2
I
Anecdote:
In 1901,Lazzerini (Italian mathematician) performed a simulation
by spinning round and dropping a needle 3407 times,He estimated
p to be 3.1415929 (accurate to the seventh number after the point!).
Exercise
Write a program to calculate p by using the method of Buffon with
d=1 and L=3/4.
Monte Carlo integration
x
d/2
L/2
qp/2
x=Lsin(q)/2
I
2/
0
s in2
p
qqdLS I
SI can be calculated from Pdp/4 (P is
the probability of finding a random
point in the area under the red curve).
Random number generation
General remark:
All the random numbers generated on a computer are in fact pseudo-
random numbers which are produced through a deterministic
algorithm!
Basic generator:
It produces uniformly distributed random numbers on [0,1],
Desirable properties of a random number generator:
It should generate sequences of random numbers which are
uniform,
uncorrelated,
with an extremely long period.
A widely used random-number generation method
Congruential method:
Xn = (aXn-1+c)mod(m) (1)
Properties of this sequence:
Xn < m
When a,c,X0 and m are integers,Eq.(1) generates integers between 0
and m-1,
Xn/m gives pseudo-random numbers distributed on (0,1).
Illustration:
With a=3,c=1,m=16 and a seed X0=2,we have the following
sequences,
2,7,6,3,10,15,14,11,2,7,…
The period of this sequence is 8.
Condition for a congruential generator to have a full period:
It can be shown that the congruential generator given in Eq.(1) has
a full period,i.e,m,if and only if
1) c and m have no common divisor (i.e.,no common divisor other
than 1);
2) a?1 mod(g) for every prime factor of m;
3) a?1 mod(4) if m is multiplier of 4.
It is obvious that m should be as large as possible! m=231-1 is often
used on a computer with 32 bits.
Statistical test of random number generators
Uniformity test,Break up the interval between 0 and 1 into a large
number of small bins and after generate a large number of random
numbers check for uniformity in the numbers of entries in each
bins (uniform histogram).
Parking lot test,Plot points in an m-dimensional space where the
m-coordinates of each point are determined by m-successive calls
to the random number generator,Then,look for regular structure
(resulting from some correlation).
Illustration:
good generator (parking lot test) bad generator(parking lot test)
1
100
1
100
How to generate distributions other than the
uniform one on (0,1)?
Transformation of random variables:
1) Uniform distribution on (a,b):
Let x be an uniform distribution on (0,1),the uniform distribution
on (a,b) is obtained from
y=a+(b-a)x
2) Random integers from 0 to N:
First divide the interval [0,1] into N+1 segments (each has a length
equal to D=1/(N+1)),M (0<M<N) is generated when
MD < x < (M+1)D
Code:
integer mran
mran=(N+1)*ran()
3) Random variables with an arbitrary distribution function
Let F(x) be the cumulative distribution function.
Remark:
F (x) is itself a random variable,Moreover,it can be shown that F(x)
is uniform random number on (0,1).
Algorithm:
1) generate a uniform random number,?1;
2)?2 = F-1(?1);
1 is a random number with the cumulative distribution F(x).
Application:
This method is most useful when F-1(x) can be obtained analytically.
Geometric illustration:
Example:
Histogram:
x
N(x)
0 1F(x) = x2
N(x)
0 1 x
N(x)
0 1
rescale the
abscissa
replot with
regular bins
1
2
x
f(x)
0
2x
xx d s sxF 2
0
2)(
Box-Müller method (generator of standard normal distribution):
Algorithm:
1) generate two uniform random numbers,?1,and?2 on (0,1);
2) n1 = (-2ln?1)1/2cos(2p?2)
n2 = (-2ln?1)1/2sin(2p?2)
give two independent random variables with standard Gaussian
distribution,i.e.,m=0 and s=1.
3) random variables with arbitrary m and s can be obtained from
n’1 = (n1 - m)/s
n’2 = (n2 - m)/s
Rejection techniques:
Basic rejection method:
Purpose:
Generate a sequence of random numbers which probability distribution
function is f(x),f(x) is bound by c and defined on a finite interval [a,b].
Algorithm:
1) generate one uniform random
number,?1,on (a,b) and one uniform
random number,?2,on (0,c);
2) accept?1 if?2? f(?1)/c;
3) repeat the above process.
Geometric meaning:
1 and?2 are the coordinates of random points over the rectangle,
When we keep only?1 for the points below f(x),?1 has the probability
distribution given by f(x).
c
a b
f(x)
x
Efficiency of rejection method:
The efficiency of a rejection is low when the surface area below f(x)
is small.
Illustration:
f1(x) can be sampled more efficiently by a rejection method that
f2(x).
c
a b
f2(x)
x
c
a b
f1(x)
x
Limitation of the method:
One must know the maximum of the distribution function and f(x)
should be defined on a finite interval.
