http://www.physics.ldu.edu.cn/kecheng/jisuanwuli/retong/my%20contents/电子教案/第七章%20Monte%20Carlo%20方法.htm
第七章 Monte Carlo 方法
稍微复杂一点的物理体系都会涉及粒子之间的相互作用，例如，一个含由N个电子组成的原子，N个原子组成的体系等。如果这些体系中每两个粒子之间的相互作用可以忽略不计，那么问题就会简单的多，我们可以对每一个粒子单独求解，这时，对这些体系的计算机处理最多只会遇到三维积分。而实际上我们需要对这些体系的更为精确的描述，这时，我们不得不面对多维积分，用我们在第一章中介绍的方法计算这些积分无疑将是非常困难的。然而Monte Carlo 方法总能非常有效地处理这类问题。
Monte Carlo方法简单地说是指用随机数来处理数值问题的一类方法。它起源于二次世界大战期间美国提出的研制原子弹的“曼哈顿计划”，该计划的主持人之一、数学家John von Neumann用摩纳哥的赌城Monte Carlo来命名这种方法。该方法的基本思想是：不用所有的求积点来计算积分的值，而只是在这些求积点中选取一些具有代表性的点，用这些点对积分值做出估计。
7.1 Monte Carlo方法计算积分
Monte Carlo方法的优势在于计算高维积分，然而现在，为了说明Monte Carlo 方法计算积分的基本思想，我们还是从较简单的一维积分说起。
我们要计算积分是
，(7.1)
将积分区间分成均匀的N等分，记
，则用矩形法计算上述积分的计算公式为：
，则用矩形法计算上述积分的计算公式为：
。(7.2)
在(7.2)式中，求积点
在求积区间[0,1]上的分布是等距的。而如果我们用Monte Carlo方法求这个积分，则是令求积点在求积区间[0,1]上随机分布。
在求积区间[0,1]上的分布是等距的。而如果我们用Monte Carlo方法求这个积分，则是令求积点在求积区间[0,1]上随机分布。
| |||
| |||
图7-1 矩形法和Monte Carlo方法求积分示意图
这样做对(7.1)式所描述的一维积分计算是没有任何好处的。因为Monte Carlo方法的计算误差由统计学规律决定，根据统计学中的中心极限定理，Monte Carlo方法的误差应为
，(7.3)
而矩形公式的误差为
。在同样的计算工作量下，矩形公式显然能够给出更加精确的结果。Monte Carlo方法的真正威力体现在高维积分的计算上，当积分维数
时，矩形公式的计算误差是
，更精确一点的Simpson公式的计算误差是
，而Monte Carlo方法的计算误差仍然是
，那么
（前者）和
（后者）就是Monte Carlo方法和另外两种计算方法的分水岭。我们在实际处理多体问题时，如果系统中的粒子数为n，则所遇到的积分的维度通常为3n，当n非常大时，Monte Carlo方法的表现将优于其他方法。
。在同样的计算工作量下，矩形公式显然能够给出更加精确的结果。Monte Carlo方法的真正威力体现在高维积分的计算上，当积分维数时，矩形公式的计算误差是，更精确一点的Simpson公式的计算误差是，而Monte Carlo方法的计算误差仍然是，那么（前者）和（后者）就是Monte Carlo方法和另外两种计算方法的分水岭。我们在实际处理多体问题时，如果系统中的粒子数为n，则所遇到的积分的维度通常为3n，当n非常大时，Monte Carlo方法的表现将优于其他方法。
下面的程序给出了用Monte Carlo方法计算一维积分的一个例子，所计算的积分为
，(7.4)
这一积分的精确解是
。在程序中，我们用Fortran内部函数RAN(I)来产生在区间[0,1)上均匀分布的随机数。我们取积分点数为N=10000，计算所得的结果为0.78461442，这一结果精确到了小数点后两位。
。在程序中，我们用Fortran内部函数RAN(I)来产生在区间[0,1)上均匀分布的随机数。我们取积分点数为N=10000，计算所得的结果为0.78461442，这一结果精确到了小数点后两位。
C SIMPLE SAMPLING MONTE CARLO
PARAMETER(N=10000)
REAL*8 SUM, X
INTEGER ISEED
SUM=0.D0
ISEED=684129
DO I=1,N
X=RAN(ISEED)
SUM=SUM+F(X)
ENDDO
SUM=SUM/N
PRINT*,N,SUM
END
DOUBLE PRECISION FUNCTION F(X)
REAL*8 X
F=1.D0/(1.D0+X*X)
END
7.1.1简单抽样的Monte Carlo方法
用Monte Carlo方法求积分，我们可能遇到两种极端情况。一种极端情况如图7-2(a)所示，被积函数为一个在求积区间的任何点取值都一样的常数，这时，我们只需要在求积区间中随意取一个点，就可以算出积分的精确值。
图7-2 用Monte Carlo方法求积分时可能遇到的两种极端情况
对于更为一般的情况，即如果被积函数在积分区间上足够光滑，则我们只需用在积分区间上均匀分布的随机数计算积分就可以了，我们称这种方法为简单抽样的Monte Carlo方法，上面程序正是用的这种方法。
7.1.2重要抽样的Monte Carlo方法
另一个极端情况如图7-2(b)所示，被积函数在某一个x0附近有一个很窄的峰，在其他x处为零，我们所选取的随机数如果在积分区间上均匀分布的话，即使选取很多的积分点，也难以将积分算准。
当遇到这种情况，我们可以考虑对Monte Carlo方法做一个改进。引入一个在积分区间内与被积函数
的行为类似的正定权函数
，其中是归一化的，
的行为类似的正定权函数，其中是归一化的，
，(7.5)
对积分做一个变形
，(7.6)
令
，(7.7)
则，
，
，
，(7.8)
因而，(7.6)式等价于
。(7.9)
(7.9)式中，被积函数变为
，而
的具体形式使得足够光滑，所以，我们可以对y在积分区间上产生均匀分布的随机数计算(7.9)式给出的积分。这时，如果我们按照(7.7)式所描述的x与y之间的关系，将在积分区间上均匀分布的y点变化成x点，那么我们就会发现，这些x点的分布是以
为概率密度的，而的行为与
相似，所以，这种选择随机数的方法实际上可以理解为在值较高的地方，选取的随机点分布的密集一些，在值较小的地方，选取的随机点分布的则比较稀疏。这样，对于，图7-2 (b)所描述的被积函数，我们可以只在
附近选取随机点，这样将大大降低达到同样计算精度所需的工作量。由于这种方法是按照它们的重要程度选取随机数，所以叫做重要抽样的Monte Carlo方法。
，而的具体形式使得足够光滑，所以，我们可以对y在积分区间上产生均匀分布的随机数计算(7.9)式给出的积分。这时，如果我们按照(7.7)式所描述的x与y之间的关系，将在积分区间上均匀分布的y点变化成x点，那么我们就会发现，这些x点的分布是以为概率密度的，而的行为与相似，所以，这种选择随机数的方法实际上可以理解为在值较高的地方，选取的随机点分布的密集一些，在值较小的地方，选取的随机点分布的则比较稀疏。这样，对于，图7-2 (b)所描述的被积函数，我们可以只在附近选取随机点，这样将大大降低达到同样计算精度所需的工作量。由于这种方法是按照它们的重要程度选取随机数，所以叫做重要抽样的Monte Carlo方法。
下面我们仍以(7.4)式的积分为例，对这种方法加以说明。选择权函数为
，由(7.7)式可以推出重要抽样的随机数x与均匀分布的随机数y之间的关系为
。下面的程序是基于这种抽样计算的积分值，计算结果为0.78532006，与简单抽样的Monte Carlo方法相比，在相同的积分点数下，将积分精度提高到了小数点后四位。
，由(7.7)式可以推出重要抽样的随机数x与均匀分布的随机数y之间的关系为。下面的程序是基于这种抽样计算的积分值，计算结果为0.78532006，与简单抽样的Monte Carlo方法相比，在相同的积分点数下，将积分精度提高到了小数点后四位。
PARAMETER(N=10000)
REAL*8 SUM, X,Y,F,W,FX,FF
INTEGER ISEED
SUM=0.D0
ISEED=684129
DO I=1,N
Y=RAN(ISEED)
X=FX(Y)
FF=F(X)/W(X)
SUM=SUM+FF
ENDDO
SUM=SUM/N
PRINT*,N,SUM
END
DOUBLE PRECISION FUNCTION F(X)
REAL*8 X
F=1.D0/(1.D0+X*X)
END
DOUBLE PRECISION FUNCTION W(X)
REAL*8 X
W=(4.D0-2.D0*X)/3.D0
END
DOUBLE PRECISION FUNCTION FX(Y)
REAL*8 Y
FX=2.d0-DSQRT(4.D0-3.D0*Y)
END
7.2 随机数的产生
从上面的描述可以看出，Monte Carlo方法计算的成功与否在很大程度上取决于所用的随机数是否随机，随机数的分布是否合适。而真正的Monte Carlo程序通常是非常耗时的，所以，如果随机数选取的不好，将会浪费掉大量的时间。
7.2.1均匀分布的伪随机数
随机数是指从某种特定分布中随机选取的数，并且每一次选取都和任何其他次的选取毫无关联。这样，当我们选了足够多的随机数后，这些随机数的分布将服从前面提到的特定分布。对于均匀分布的随机数，每个数出现的概率都与任何其他数出现的概率相等。
计算机上产生的随机数都是基于某种算法，也就是说，什么时候出现什么样的数是完全可预测的，但是，对于不知道这种算法的人而言，这些数看起来是随机的，我们称它为“伪随机数”。
产生均匀分布的伪随机数是一种基本的计算机操作，它是产生其它分布的随机数的基础。一种最常用的产生均匀分布的伪随机数的计算机算法是线性同余法：
，(7.10)
式中的a、b、M是正整数称为幻数，mod意味着取余。第一个数是人为给出的，称为种子，从一个种子出发，可以得到一个伪随机数序列。例如，令a=6，b=7，M=5，产生的序列为{3，0，2，4，1，3，…}；令a=7，b=0，M=10，产生的序列为，{3，1，7，9，3，…}。
从这两个例子可以看出，线性同余法产生的伪随机数序列分布于[0，M-1]之间，并以不大于M的周期重复。前面的两个例子的重复周期都很短，显然，所产生的序列不够“随机”。
通常，对于计算机产生的为随机数序列，我们有如下要求：
1， 重复周期长；
2， 数之间的关联小；
3， 计算速度快。
对于线性同余法，通常选择a=16807，b=0，M=231-1，已经验证，这时产生的随机数序列的重复周期为M=231-1，数的分布足够随机，并且算法速度很快。如果要产生[0,1)之间的随机数，则可以令所产生的随机数序列中的每一个数除以M。
我们可以用Fortran内部函数RAN(ISEED)、内部子程序RANDOM_NUMBER、或运行时间子程序RANDOM来产生在区间[0,1)上均匀分布伪随机数。所有的函数、子程序所用的都是基于线性同余法的同一种算法，所以，只要给出同一个种子，它们返回的随机数序列将是一样的。下面是用调用内部子程序RANDOM_NUMBER来产生10个随机数组成的序列的一个例子。程序中“CALL SEED(7531)”语句产生随机数序列的种子，如果要求每次运行时返回不同的随机数，只需将其替换为“CALL SEED(RND$TIMESEED)”。
USE DFLIB
REAL*8 RAND
OPEN(FILE='RAND.DAT',UNIT=10)
CALL SEED(7531)
DO I=1,10000
CALL RANDOM_NUMBER (RAND)
WRITE(10,*)I,RAND
ENDDO
END
7.2.2特定分布的伪随机数
有了均匀分布的随机数，我们就可以用它们产生具有某种特定分布的随机数，例如，Poisson分布
，(7.11)
或Gaussian分布
。(7.12)
假定
为在区间[x,x+dx]找到随机变量x的概率，做一个变量替换
，在区间[y,y+dy]找到变量y的概率记为p(y)dy，则有
为在区间[x,x+dx]找到随机变量x的概率，做一个变量替换，在区间[y,y+dy]找到变量y的概率记为p(y)dy，则有
。(7.13)
如果y在区间[0,1)上均匀分布，
，(7.14)
则有
。（7.15）
如果为Poisson分布，则有
为Poisson分布，则有
。（7.16）
令，由(7.16)式求反函数，得
，由(7.16)式求反函数，得
。(7.17)+
如果
是在[0,1)上均匀分布的随机数，则
的分布则满足Poisson概率分布函数(7.11)。下面的程序利用计算机内部函数产生在[0,1)区间上均匀分布的随机数，然后用上述算法，产生在
上按Poisson分布函数分布的随机数。所得的结果如图7-3所示。
是在[0,1)上均匀分布的随机数，则的分布则满足Poisson概率分布函数(7.11)。下面的程序利用计算机内部函数产生在[0,1)区间上均匀分布的随机数，然后用上述算法，产生在上按Poisson分布函数分布的随机数。所得的结果如图7-3所示。
图7-3 由文中给出的程序计算得到的(a)在区间[0,1)上均匀分布的随机数；(b)在区间上按Poisson分布函数分布的随机数。
上按Poisson分布函数分布的随机数。
PARAMETER(N=10000)
REAL*8 X,Y,FX
INTEGER ISEED
ISEED=684129
OPEN(FILE='RAND.DAT',UNIT=10)
DO I=1,N
Y=RAN(ISEED)
X=FX(Y)
WRITE(10,*)I,Y,X
ENDDO
END
DOUBLE PRECISION FUNCTION FX(Y)
REAL*8 Y
FX=-LOG(1.D0-Y)
END
7.2.3 Metropolis算法
在上面的例子中，为了产生具有特定分布的随机数序列，我们需要知道的解析形式，对其做积分，并求出反函数
。为了将被积函数的行为描述的更好，的形式可能会非常复杂，这时，虽然我们可以用数值方法实现上面的步骤，但这会带来不必要的麻烦。
的随机数序列，我们需要知道的解析形式，对其做积分，并求出反函数。为了将被积函数的行为描述的更好，的形式可能会非常复杂，这时，虽然我们可以用数值方法实现上面的步骤，但这会带来不必要的麻烦。
Metropolis算法则可以很方便地产生具有任意形状概率分布的随机数。下面我们不加证明地对Metropolis算法进行描述。我们不再将考虑的范围局限于一维情况，而是考虑n维位形空间中的概率密度函数
，其中
为n维空间中的矢量。
，其中为n维空间中的矢量。
假设一个随机行走者在n维空间中运动，它从空间中的某一点
开始，相继走出各步的终点产生出一个序列
，如果走的步数足够多，那么这些点列的分布将满足概率密度函数。
开始，相继走出各步的终点产生出一个序列，如果走的步数足够多，那么这些点列的分布将满足概率密度函数。
随机行走者遵从的规则如下：假设行走者处于
上，它将试探性地迈出一步到达一个新点
上，它将试探性地迈出一步到达一个新点
，(7.18)
这个新点可以用任何方便的方法选取，例如，可以在周围边长为h的小立方体中均匀地选取。这时，可令
的分量
周围边长为h的小立方体中均匀地选取。这时，可令的分量
，(7.19)
其中，
为在[0,1]上均匀分布的随机数。新点
是否被接受将取决于比值
为在[0,1]上均匀分布的随机数。新点是否被接受将取决于比值
，(7.20)
如果
，那么接受这一步，即令
；如果
，那么就以概率
接受这一步，即将与在区间[0,1]上均匀分布的随机数
比较，若
，就接受这一步，否则舍弃这一步，舍弃这一步意味着
。然后按照同样的规则产生
，
，……。
，那么接受这一步，即令；如果，那么就以概率接受这一步，即将与在区间[0,1]上均匀分布的随机数比较，若，就接受这一步，否则舍弃这一步，舍弃这一步意味着。然后按照同样的规则产生，，……。
Metropolis算法在实际应用时一个需要注意的问题就是步长h的选取就。我们注意到，算法所产生各点彼此间不是相互独立的，
的产生依赖于点的位置。假设是使取极大值的一个点，如果h取得过大，那么
可能比
小得多，大部分的试验步将不被接受，因此抽样的效率很低；如果h取得过小，以后产生的诸点将分布在点附近，点列的分布与所的描述符合得不好；通常的办法是令h的选择使得试验步被接受的几率为
。
的产生依赖于点的位置。假设是使取极大值的一个点，如果h取得过大，那么可能比小得多，大部分的试验步将不被接受，因此抽样的效率很低；如果h取得过小，以后产生的诸点将分布在点附近，点列的分布与所的描述符合得不好；通常的办法是令h的选择使得试验步被接受的几率为。
下面的程序将应用Metropolis算法产生具有Poisson分布的伪随机数。
PARAMETER(N=10000)
REAL*8 X,W,H,RND
H=0.5D0
X=0.7D0
OPEN(FILE='RAND.DAT',UNIT=10)
DO I=1,N
CALL RANDOM_NUMBER(RND)
CALL METROPOLIS(X,H,RND)
WRITE(10,*)I,X
ENDDO
END
SUBROUTINE METROPOLIS(X,H,RND)
REAL*8 X,H,XT,RND
XT=X+H*(2.D0*RND-1.D0)
IF(XT.LT.0.D0) RETURN
R=W(XT)/W(X)
IF(RND.GT.R)RETURN
X=XT
END
DOUBLE PRECISION FUNCTION W(X)
REAL*8 X
W=EXP(-X)
END
No comments:
Post a Comment