Monte Carlo Methods

1
Monte Carlo Methods

• So far we have discussed Monte Carlo methods
  based on a uniform distribution of random numbers
  on the interval 0,1
• p(x) 1 0 ? x ? 1 p(x) 0
  otherwise
• the hit and miss algorithm generates pairs of
  points (x,y) and either accepts or rejects the
  point
• the sample mean method samples points uniformly
  on the interval a,b

2
(No Transcript)
3
(No Transcript)
4
Choose p(x) so that f(x)/p(x) is a fairly
constant function
5
Nonuniform probability distributions

• Consider a probability density p(x) such that
  p(x)dx is the probability that event x is in the
  interval between x and xdx
• ? We have
  ? p(x) dx 1
  -?
• x Calculate P(x)
  ? p(x) dx r ? dr
  -?
• p(x)dx dr gt p(x) dr/dx
• r is a uniform random number on the interval
  0,1
• Invert this and solve for x in terms of r

6
Eg.1 suppose we want p(x) 1/(b-a), a? x ? b
0
otherwise

• x
• Calculate P(x) ? p(x) dx
• a
• r (x-a)/(b-a)
• Solving for x, x a (b-a) r
• obvious!

7
Eg.2 p(x) (1/?) exp(-x/?), 0,? 0
x lt 0

• Hence r P(x) 1 - exp(-x/?)
• And x -? ln(1-r) -? ln(r)
• This technique is only feasible if the inversion
  process can be carried out.

8
Eg.3 p(x) (1/ 2??2)1/2 exp(-x2/2?2)
P(x) ?

• However the two-dimensional distribution
• p(x,y)dxdy (1/ 2??2) exp(-(x2y2)/2?2)dxdy can
  be integrated by a change of variable
• ? r2/2 (x2y2)/2?2 , tan?
  y/x
• p(?,?)d?d? (1/2?) exp(-?)d?d?
• Generate ? uniformly on the interval 0,2?
  i.e. ? 2? r
• Generate ? according to the exponential
  distribution with ? 1 i.e. ? -ln r
• x (2??2)1/2 cos? and y(2??2)1/2 sin? are
  Gaussian distributed

9
Importance Sampling

• The error estimate in Monte Carlo is proportional
  to the variance of the integrand
• (ltf2gt - ltfgt2 )1/2
  
  n1/2
• How can we reduce the variance?
• 
  b Introduce a function p(x) such that
  ? p(x)dx 1
  
  a
• b
  And rewrite F ? f(x)/p(x) p(x) dx
  a

10
Importance Sampling

• Evaluate the integral by sampling according to
  p(x)
• Fn (1/n) ? f(xi )/p(xi )
• Choose p(x) to minimize variance of f(x)/p(x)
• i.e try to make f(x)/p(x) slowly varying since a
  constant has zero variance

11
1eg. F ? exp(-x2) dx
.746824 0
Sample 0,1 uniformly n Fn
?n 100.
0.74613 0.19602 1000.
0.75169 0.19673 10000. 0.74812
0.19993 Choose p(x) A exp(-x) and sample
0,1 again n Fn
?n 100. 0.75105
0.05076 1000. 0.74882 0.05411
10000. 0.74753 0.05420 The
variance of the integrand is reduced by about a
factor of 4 which means that fewer samplings are
needed to obtain the same accuracy.
12
c Program Importance Sampling n10000
h1. a0. b1. sum0.
psum0. sum20. psum20. m2
do 4 i1,n wwr250(idum) xabww
y-log(1.-wwwwexp(-1.)) gexp(-yy)
p(1.-wwwwexp(-1.))/(1.-exp(-1.))
fexp(-xx) sumsumf psumpsumg/p
sum2sum2ff psum2psum2(gg)/(pp)
sig2sum2/i -(sum/i)(sum/i) sigsqrt(sig2)
psig2psum2/i -(psum/i)(psum/i)
psigsqrt(psig2) if((i-(i/10m)10m).eq.0)
then write(6,10) 1.i,sum/i,sig,psum/i,psig
mm1 else continue end if 10
format(1x,f10.0,3x,f10.5,3x,f10.5,3x,f10.5,3x,f10.
5) 4 continue stop end
13
The above ideas can be used to simulate many
different types of physical problems

• Random walks and polymers
• Percolation
• Fractal growth
• Complexity and neural networks
• Phase Transitions and Critical Phenomena

14
Monte Carlo Methods

• Random numbers generated by the computer are used
  to simulate naturally random processes
• many previously intractable thermodynamic and
  quantum mechanics problems have been solved using
  Monte Carlo techniques
• how do we know is the random numbers are really
  random?

15
Random Sequences

• A sequence of numbers r1,r2, is random if there
  are no correlations among the numbers in the
  sequence
• however most random number generators yield a
  sequence in which each number is used to find the
  succeeding one according to a well defined
  algorithm
• the most widely used random number generator is
  based on the linear congruential method

16
Given a seed x0, each number in the sequence is
determined by the one-dimensional map

• where a,c and m are integers
• the notation y z mod m means that m is
  subtracted from z until 0? y ltm
• the process is characterized by the multiplier a,
  the increment c and the modulus m
• since m is largest integer generated by this
  method, the maximum possible period is m

17
Example

• a3 c4 m32 and x01 produces
• x1(3 x 1 4) mod32 7
• x2(3 x 7 4) mod32 25
• x3(3 x 25 4) mod32 79mod3215
• and so on .
• 1,7,25,15,17,23,9,31,1,7,25,.
• period is 8!
• Rather than the maximum of 32

18
Random Sequences

• If we choose a, c and m carefully then all
  numbers in the range from 0 to m-1 will appear in
  the sequence
• to have the numbers in the range 0 ? r lt1, the
  generator returns xm/m which is always lt 1
• there is no necessary and sufficient test for the
  randomness of a finite sequence of numbers
• we need to consider various tests
• an obvious requirement for a random number
  generator is that its period be much greater than
  the number of random numbers needed in a specific
  problem

19
Sequences

• A way of visualizing the period is to consider a
  random walker and plot the displacement as a
  function of the number of steps N
• when the period of the random number generator is
  reached the plot will begin to repeat itself
• consider a899, c0, m32768 with x012

Sequence
20
Correlations

• We can check for correlations by plotting xik as
  a function of xi
• if there are any obvious patterns in the plot
  then there are correlations

correlations
21
Uniformity
test
22
Correlations

• One way to reduce sequential correlation and to
  lengthen the period is to mix or shuffle two
  different random number generators
• statistical tests should be performed on random
  number generators for serious calculations

23
Exploration
xi1 4xi(1 - xi)

• It has been claimed that the logistic map in the
  chaotic region is a good random number generator
• test this for yourself

• will make the sequence uniform