Monte Carlo Method (2)
Monte Carlo Integration
Hit and miss method:
b
a
xd xgI )(
When -d? g(x)?c,we can always first make a shift f(x)=g(x)+d
so that f(x)? 0,
Algorithm:
1) generate N uniform random numbers,?i (i=1,2,…,N),on (a,b)
and N uniform random numbers,?i (i=1,2,…,N),on (0,c);
2) calculate g(?i) and count the number of cases,Nhit,for which
i < g(?i);
3) I can be estimated by
Iest=c(b-a)Nhit/N
0? g(x)? c
Geometrical illustration:
c
a b x
g(x)
g(x)
miss
hit
)(
)(
abc
xd x g
p
b
ah it
N
N
Variance:
IabcNIpabcI e s t )()?v a r ()()v a r ( 2
Remark:
Nhit has a binomial distribution.
var(Iest) increases with c(b-a)-I (i.e.,the numerical error of the
method increases with c(b-a)-I).
Simple Mean-Value Method:
Remark:
<g(x)> is the mean value of g(x) for a uniform random distribution
on (a,b).
Algorithm:
generate N random numbers,xi,distributed uniformly on (a,b);
I can be estimated by using,
)()(1)()()( xgababxd x gabxd x gI
b
a
b
a
Ni ie s t xI gN ab 1 )(
Variance:
))(va r ())(va r ()()(va r)va r ( )()(
2
12
2
1 xN
abx
N
abxI ggg
N
ab N
i i
N
i ie s t
b
a
Ig xdxabN 22 )()(1
Efficiency of Monte Carlo Integration
Smaller is the variance of the estimated result,higher is the efficiency
of the method.
Illustration:
var(Iest)~1/N
For fixed N,smaller variance leads to smaller calculation error.
N
Iest
I
N
Iest
I
large variance small variance
Comparison of the efficiency of hit-miss method
with that of simple mean-value method
1) Hit-miss method
II
b
am i s sh i t
xdx gabcNe s t 2)()(1)v a r(
2) Simple mean-value method:
II
b
as i m p l e
xgxdx gabNe s t 2)()()(1)v a r(
Since c?g(x),var(Iest)hit-miss> var(Iest)simple,The simple mean-value
method is more efficient than the hit-miss method.
Remark:
Different Monte Carlo integration methods do not have the same
efficiency.
How to improve the efficient of Monte Carlo
integration method?
One widely used method is based on the variance reduction.
Importance sampling method:
Basic ideal:
fxf
xgxf
xf
xgdxxd xgI
)(
)()(
)(
)()( f(x) > 0
Intuitively,one expects a small variance when g(x)/f(x) is nearly
constant,
Variance:
Ig xf xdxxf xg 2
2
)(
)(
)(
)(v a r
When g(x)>0 and
)()()( xdx gxgxf
,var[g(x)/f(x)]=0.
In more general cases,when
)(
)()(
xgdx
xgxf,var[g(x)/f(x)] reaches a
minimum.
Geometrical implication:
Simple Mean-Value Method Important Sampling Method
x
g(x)
0
<g(x)>
x
g(x)/f(x)
0
f(x)?g(x)
f(x)
)( )(xf xg
MR2T2 Algorithm (Metropolis Algorithm):
Reference:
N,Metropolis,A.W,Rosenbluth,M.N,Rosenbluth,A.H,Teller and
E,Teller,J,Chem,Phys,21,1087,(1953).
Typical integrals often encountered in statistical mechanics:
Example,Internal energy
<E (rN)> =?drN E(rN)r(rN)
r(rN)=exp(- E(rN)/kT)/Q
Q =?drN exp(-E(rN)/kT)
Remark:
r(rN) is distribution function with a narrow peak near <E (rN)>.
This is why a simple mean-value method will fail badly for
calculating accurately <E (rN)> and Q.
Basic idea:
Generate more points in the region where the probability density
function,pdf,is large.
Difficulty:
Since we do not know completely the distribution function (unknown
normalization constant),it is impossible to use the usual techniques,
e.g.,transformation method,rejection method.
Solution:
Metropolis et al proposed an iterative method to circumvent the
problem,
Sketch of MR2T2 algorthm for sampling a random variable with an
arbitrary pdf,f(x):
1) Start with random number,?0,which can have any distribution
function;
2) Make an attempt for generating a next random number by
’i=?i-1 +h and h is an uniform random number on (-d,d);
3) If f(?’i)>f(?i-1),?’1 is accepted (i.e.,?i=?’i),If f(?’i)<f(?i-1),is
accepted with a probability equal to f(?’i)/f(?i-1) when?’i is not
accepted,?i=?i-1;
4) Go to (2) and repeat the process.
The distribution of the random numbers generated by this algorithm
tends to be f(x) asymptotically.
Caution:
1) When f(x) is defined on a finite interval,e.g.,(a,b),if the
generated?’i is outside (a,b),step (2) should be continued until a
’i inside (a,b) is obtained,It is incorrect to merely keep?i-1 (i.e.,
i=?i-1) when?’i is outside (a,b).
2) In step (2),?’i must be generated from a symmetric distribution
centered at?i-1,
Influence of the parameter d:
When d is small,the acceptation ratio is high but the successively
generated random numbers evolves slowly,When d is too large,
the acceptation ratio will become very low and the algorithm loses
the efficient.
Practical recipe for the value of d:
Adjust d so that the acceptation ration is equal to 0.5.
Illustration with a simple example:
Generate a sequence of random numbers which take only the values
1 and 2 with f(2)/f(1)=3.
The following sequence is generated with?0=1,h=+1 or -1 and
p+=p-=0.5.
122222122122122121222221222221 8 22
212121212222121222121212222122 10 20
122222222212212221221212121222 8 22
221212122222222212122121222222 7 23
221222221222222222221212122222 5 25
38 112 112/38=2.94
MR2T2 algorithm as a stochastic process
Many formal properties of algorithm can be proven by using the tools
of stochastic processes.
Basic idea:
Consider the successive steps of a random number or more generally
a random configuration generation as a stochastic process,
A stochastic process describes the time evolution of a system with
probabilities.
Main quantities for characterizing a stochastic process:
Consider a process with N states (i=1,2,…N).
pi(t),probability for the system being in the state i at t.
p(j,t+1|i,t;kt-1,t-1;…;k 0,0),transition probability.
Markovian process (Markov chain):
p(j,t+1|i,t;kt-1,t-1;…;k 0,0)=p(j,t+1|i,t)
p(j,t+1|i,t) is also often simply denoted as pij.
Master equation:
The time evolution of a Markovian process is described by the following
master equation,
In the investigation of stochastic processes,one knows the transition
probability and studies the time evolution of pi(t).
When devising a sampling method,one looks for a suitable transition
probability so that pi(t) tends asymptotically to a specified distribution.
The stationary solution of Eq.(1) is given by
ij ijiij jiji pp ttdt td )(_)()( ppp (1)
(2)
ij ijiij jij pp
0_ pp
Detailed balance:
pipij = pjpji
Other conditions to be satisfied by pij:
pij > 0
j ijp 1
j jijpi pp
The transition probability proposed for a sampling algorithm must
satisfy these conditions and the detailed balance.
MR2T2 algorithm in terms of transition probability
p*ij pj?pi,j?i
p*ij pj/pi pj<pi,j?i
Pij =
ij ijii pp 1
Usually,the transition is made in two step,p*ij is the probability of
the attempt transition.
Important note:
p*ij must be symmetric,i.e.,p*ij = p*ji to ensure the process is ergotic
(i.e.,all the states can be reached).
It is easy to check that the transition probability prescribed above
satisfies the detailed balance and all the other conditions,
(pii,probability of staying in the same state)
How a system evolves toward the prescribed
distribution in MR2T2 algorithm?
Consider a very large ensemble of systems,At time t,there are ni(t)
systems in the state i and nj(t) systems in the state j,Suppose pj> pi,
The number of systems moving from j to i is nj(t)p*ji and that from i
to j is ni(t)p*ijpi/ pj.
The net number of systems moving from i to j is
ni?j= p*ij[ni(t)pi/pj-nj(t)] (p*ji= p*ij is used)
If ni(t) /nj(t) > pi/pj,?ni?j>0.
On average,there are more systems move from i to j,This makes
the distribution of the systems in the ensemble tends to the
prescribed distribution,p,
Advantages and drawbacks of MR2T2 algorithm
Advantage:
1) robust and versatile (can be used for any distribution even with
unknown normalization constant);
2) simple to implement.
Drawbacks:
1) iterative procedure;
2) there is some correlation between the successive configurations (or
random numbers) generated by this method.
Remark:
When possible,it is preferable to use a method which generates
independent random numbers.
Exercise
Write a program to generate a standard Gaussian distribution
from an uniform distribution on (0,1),Adjust d so that the
acceptation ratio is 0.5,Discard the random numbers
generated in the initial transition period and see how long is
the transition period.
Monte Carlo Integration
Hit and miss method:
b
a
xd xgI )(
When -d? g(x)?c,we can always first make a shift f(x)=g(x)+d
so that f(x)? 0,
Algorithm:
1) generate N uniform random numbers,?i (i=1,2,…,N),on (a,b)
and N uniform random numbers,?i (i=1,2,…,N),on (0,c);
2) calculate g(?i) and count the number of cases,Nhit,for which
i < g(?i);
3) I can be estimated by
Iest=c(b-a)Nhit/N
0? g(x)? c
Geometrical illustration:
c
a b x
g(x)
g(x)
miss
hit
)(
)(
abc
xd x g
p
b
ah it
N
N
Variance:
IabcNIpabcI e s t )()?v a r ()()v a r ( 2
Remark:
Nhit has a binomial distribution.
var(Iest) increases with c(b-a)-I (i.e.,the numerical error of the
method increases with c(b-a)-I).
Simple Mean-Value Method:
Remark:
<g(x)> is the mean value of g(x) for a uniform random distribution
on (a,b).
Algorithm:
generate N random numbers,xi,distributed uniformly on (a,b);
I can be estimated by using,
)()(1)()()( xgababxd x gabxd x gI
b
a
b
a
Ni ie s t xI gN ab 1 )(
Variance:
))(va r ())(va r ()()(va r)va r ( )()(
2
12
2
1 xN
abx
N
abxI ggg
N
ab N
i i
N
i ie s t
b
a
Ig xdxabN 22 )()(1
Efficiency of Monte Carlo Integration
Smaller is the variance of the estimated result,higher is the efficiency
of the method.
Illustration:
var(Iest)~1/N
For fixed N,smaller variance leads to smaller calculation error.
N
Iest
I
N
Iest
I
large variance small variance
Comparison of the efficiency of hit-miss method
with that of simple mean-value method
1) Hit-miss method
II
b
am i s sh i t
xdx gabcNe s t 2)()(1)v a r(
2) Simple mean-value method:
II
b
as i m p l e
xgxdx gabNe s t 2)()()(1)v a r(
Since c?g(x),var(Iest)hit-miss> var(Iest)simple,The simple mean-value
method is more efficient than the hit-miss method.
Remark:
Different Monte Carlo integration methods do not have the same
efficiency.
How to improve the efficient of Monte Carlo
integration method?
One widely used method is based on the variance reduction.
Importance sampling method:
Basic ideal:
fxf
xgxf
xf
xgdxxd xgI
)(
)()(
)(
)()( f(x) > 0
Intuitively,one expects a small variance when g(x)/f(x) is nearly
constant,
Variance:
Ig xf xdxxf xg 2
2
)(
)(
)(
)(v a r
When g(x)>0 and
)()()( xdx gxgxf
,var[g(x)/f(x)]=0.
In more general cases,when
)(
)()(
xgdx
xgxf,var[g(x)/f(x)] reaches a
minimum.
Geometrical implication:
Simple Mean-Value Method Important Sampling Method
x
g(x)
0
<g(x)>
x
g(x)/f(x)
0
f(x)?g(x)
f(x)
)( )(xf xg
MR2T2 Algorithm (Metropolis Algorithm):
Reference:
N,Metropolis,A.W,Rosenbluth,M.N,Rosenbluth,A.H,Teller and
E,Teller,J,Chem,Phys,21,1087,(1953).
Typical integrals often encountered in statistical mechanics:
Example,Internal energy
<E (rN)> =?drN E(rN)r(rN)
r(rN)=exp(- E(rN)/kT)/Q
Q =?drN exp(-E(rN)/kT)
Remark:
r(rN) is distribution function with a narrow peak near <E (rN)>.
This is why a simple mean-value method will fail badly for
calculating accurately <E (rN)> and Q.
Basic idea:
Generate more points in the region where the probability density
function,pdf,is large.
Difficulty:
Since we do not know completely the distribution function (unknown
normalization constant),it is impossible to use the usual techniques,
e.g.,transformation method,rejection method.
Solution:
Metropolis et al proposed an iterative method to circumvent the
problem,
Sketch of MR2T2 algorthm for sampling a random variable with an
arbitrary pdf,f(x):
1) Start with random number,?0,which can have any distribution
function;
2) Make an attempt for generating a next random number by
’i=?i-1 +h and h is an uniform random number on (-d,d);
3) If f(?’i)>f(?i-1),?’1 is accepted (i.e.,?i=?’i),If f(?’i)<f(?i-1),is
accepted with a probability equal to f(?’i)/f(?i-1) when?’i is not
accepted,?i=?i-1;
4) Go to (2) and repeat the process.
The distribution of the random numbers generated by this algorithm
tends to be f(x) asymptotically.
Caution:
1) When f(x) is defined on a finite interval,e.g.,(a,b),if the
generated?’i is outside (a,b),step (2) should be continued until a
’i inside (a,b) is obtained,It is incorrect to merely keep?i-1 (i.e.,
i=?i-1) when?’i is outside (a,b).
2) In step (2),?’i must be generated from a symmetric distribution
centered at?i-1,
Influence of the parameter d:
When d is small,the acceptation ratio is high but the successively
generated random numbers evolves slowly,When d is too large,
the acceptation ratio will become very low and the algorithm loses
the efficient.
Practical recipe for the value of d:
Adjust d so that the acceptation ration is equal to 0.5.
Illustration with a simple example:
Generate a sequence of random numbers which take only the values
1 and 2 with f(2)/f(1)=3.
The following sequence is generated with?0=1,h=+1 or -1 and
p+=p-=0.5.
122222122122122121222221222221 8 22
212121212222121222121212222122 10 20
122222222212212221221212121222 8 22
221212122222222212122121222222 7 23
221222221222222222221212122222 5 25
38 112 112/38=2.94
MR2T2 algorithm as a stochastic process
Many formal properties of algorithm can be proven by using the tools
of stochastic processes.
Basic idea:
Consider the successive steps of a random number or more generally
a random configuration generation as a stochastic process,
A stochastic process describes the time evolution of a system with
probabilities.
Main quantities for characterizing a stochastic process:
Consider a process with N states (i=1,2,…N).
pi(t),probability for the system being in the state i at t.
p(j,t+1|i,t;kt-1,t-1;…;k 0,0),transition probability.
Markovian process (Markov chain):
p(j,t+1|i,t;kt-1,t-1;…;k 0,0)=p(j,t+1|i,t)
p(j,t+1|i,t) is also often simply denoted as pij.
Master equation:
The time evolution of a Markovian process is described by the following
master equation,
In the investigation of stochastic processes,one knows the transition
probability and studies the time evolution of pi(t).
When devising a sampling method,one looks for a suitable transition
probability so that pi(t) tends asymptotically to a specified distribution.
The stationary solution of Eq.(1) is given by
ij ijiij jiji pp ttdt td )(_)()( ppp (1)
(2)
ij ijiij jij pp
0_ pp
Detailed balance:
pipij = pjpji
Other conditions to be satisfied by pij:
pij > 0
j ijp 1
j jijpi pp
The transition probability proposed for a sampling algorithm must
satisfy these conditions and the detailed balance.
MR2T2 algorithm in terms of transition probability
p*ij pj?pi,j?i
p*ij pj/pi pj<pi,j?i
Pij =
ij ijii pp 1
Usually,the transition is made in two step,p*ij is the probability of
the attempt transition.
Important note:
p*ij must be symmetric,i.e.,p*ij = p*ji to ensure the process is ergotic
(i.e.,all the states can be reached).
It is easy to check that the transition probability prescribed above
satisfies the detailed balance and all the other conditions,
(pii,probability of staying in the same state)
How a system evolves toward the prescribed
distribution in MR2T2 algorithm?
Consider a very large ensemble of systems,At time t,there are ni(t)
systems in the state i and nj(t) systems in the state j,Suppose pj> pi,
The number of systems moving from j to i is nj(t)p*ji and that from i
to j is ni(t)p*ijpi/ pj.
The net number of systems moving from i to j is
ni?j= p*ij[ni(t)pi/pj-nj(t)] (p*ji= p*ij is used)
If ni(t) /nj(t) > pi/pj,?ni?j>0.
On average,there are more systems move from i to j,This makes
the distribution of the systems in the ensemble tends to the
prescribed distribution,p,
Advantages and drawbacks of MR2T2 algorithm
Advantage:
1) robust and versatile (can be used for any distribution even with
unknown normalization constant);
2) simple to implement.
Drawbacks:
1) iterative procedure;
2) there is some correlation between the successive configurations (or
random numbers) generated by this method.
Remark:
When possible,it is preferable to use a method which generates
independent random numbers.
Exercise
Write a program to generate a standard Gaussian distribution
from an uniform distribution on (0,1),Adjust d so that the
acceptation ratio is 0.5,Discard the random numbers
generated in the initial transition period and see how long is
the transition period.
