随机模拟与蒙特卡罗法

流星赶月0454

分子聚集体或宏观物系的结构与物理化学性质可以通过求相应物理量的系综平均值得
到。但是，系统的
状态数常常是一个近乎天文的数字。比如，对一个大小的一个二维点阵，即使每个点只能
取两个状态，
系统可能取的状态总数为。这样，求一个平均值，就是用速度为每秒数万亿次的巨型计算
机，
它也要算上至少年。显然，这是根本无法做到的。计算机随机模拟方法是用一小部分“代
表”状态上
的算术平均值近似系综平均值。
    计算机随机模拟是在电子计算机上对随机现象进行模拟并进而得到问题解答的方法的
简称。其中所说的
随机现象可以是原问题本身所固有的，也可以是人为建立起来的与原问题答案有一定关系
的某随机现象。计
算机随机模拟方法亦称蒙特卡罗方法(Monte Carlo Method)，后一名称已被国内外广大学
者所普遍采用。蒙
特卡罗方法的程序大都为专门应用设计，一般的应用程序不会太复杂，可以自己编写。一
个模拟退火算法程
序ASA(Adaptive Simulated Annealing)
可以在网上免费得到(http://www.ingber.com，
ftp://ftp.ingber.com)，用于许多优化的目的。
1. 蒙特卡罗方法的名称与历史
    蒙特卡罗方法为计算数学中的一种计算方法，它的基本特点是，以概率与统计中的理
论与方法为基础，
以是否适于在计算机上使用为重要标志。因此，它虽属计算方法但又与一般计算方法有很
大区别。蒙特卡罗
是摩纳哥的一个著名城市，以赌博闻名于世。蒙特卡罗方法借用这一城市的名称，属于象
征性的，是为了表
明该方法的上述基本特点的。蒙特长罗方法也称统计试验方法或计算机随机模拟方法，这
些名称同样是为表
明该方法的上述基本特点的。   
    蒙特卡罗方法作为一种可行的计算方法，是由Ulam和Von Neumann在20世纪40年代中
叶为解决研制核武器
中的计算问题而首先提出并加以运用的。在此之前，作为该方法的基本思想，实际上早已
被统计学家所发现
和利用了。例如，早在17世纪的时候，人们就知道了依频数来决定概率的方法。又如，在
19世纪末，曾有很
多人进行随机投针试验。根据针与地面上平行线束（距离均为二倍针长）相交的概率P等
于的倒数，用频数
n/ N替代概率P，并进而得到。为了使p的有效数字达到4位，置信水平为0.95，所需投针
次数要
在40万以上。因此，在还不具备实现这样大量试验的条件之前，除非为其他目的，如上例
求p是为了验证大
数定律，不会有人用进行实际试验的办法来计算所要计算的值。
    进入20世纪40年代中叶，出现了电子计算机，使得用数学方法在电子计算机上模拟这
样大量的试验成为
可能。另外，科学技术的不断发展，出现了越来越多的复杂而困难的问题，用通常的解析
方法或数值方法都
很难得到解决。蒙特卡罗方法就是在这些情况下，作为一种可行的而且是不可缺少的计算
方法被提出和迅速
发展起来的。

2. 蒙特卡罗方法的一般原理
    蒙特卡罗方法解题的一般过程是，首先构成一个概率空间；然后在该概率空间中确定
一个依赖随机变量
A(可以为任意维)的统计量g(x)，其数学期望

    作为G的近似估计。
    由以上过程看出，蒙特卡罗方法解题的最关键一步是，确定一个统计量，其数学期望
正好等于所要求的
值。为方便起见，后面将总称这样的统计量为无偏统计量。
    由于其他原因，如确定数学期望为G的统计量g(x)有困难，或为其他目的，蒙特卡罗
方法有时也用G的渐
近无偏估计代替一般过程中的无偏估计，并用此渐近无偏估计作为G的近似估计。
    蒙特卡罗方法的最低要求是，能确定这样一个与计算步数N有关的统计估计量
，当时,依
概率收敛于所要求的值G，亦即，对于任意的＞0，应有
   
    其中P(·)表示事件·的概率。
子样：具有已知分布f(x)的总体，作为它的子样形式，在蒙特卡罗方法中通常被大家所采
用的是简单子样。
简单子样指的这样的随机子样，它们之间相互独立，对于任意的x和所有的n = 1，…，
N满足
   
    简单子样具有三个非常重要的性质。第一个性质是，对于任意的统计量具有同一分布
且
相互独立，数学期望均等于

蒙特卡罗方法的收敛性：蒙特卡罗方法的近似估计依概率1收敛于G，亦即
   
    的充分必要条件是无偏统计量g(x)满足
   
    就是说，蒙特卡罗方法的收敛性取决于所确定的无偏统计量是否绝对可积，确切些说
，是它的绝对数学
期望是否存在。
    如果无偏统计量g(x)满足条件

    亦即依概率1收敛G的速度为。就是说，蒙特卡罗方法出收敛速度取决于所确定的无偏
统计量是
几次绝对可积，但收敛速度总不会超过。
蒙特卡罗方法的误差：根据中心极限定理，只要所确定的无偏统计量具有有限的异于零的
方差，对于任
意非负的X均有
   
因此，当N足够大时，就可以认为有如下近似等式：
   
    其中为置信度，1 - 为置信水平。
    根据以上结果，我们可以根据问题的要求，确定出置信水平，然后按照正态积分表确
定出X来，而近似
估计 与真值G之间的误差就可以用下式近似地得到：   
      
    通常取X为0.6745、1.96或3，相应的置信水平依次为0.5、0.95或0.997。
    蒙特卡罗方法误差容易确定。对于－般计算方法，要想给出计算结果与真值的误差并
不是一件容易办到
的事情，蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，其中的X是确定的，N是实
际的抽样数，未知
的仅仅是均方差。可是，它是可以通过计算的同时计算给出的，即
      
    对于再复杂的蒙特卡罗方法计算问题，均方差的确定大致上也是如此，因此，误差容
易确定成了蒙特卡
罗方法另一大优点。
    一般计算方法常存在着有效位数损失问题，要想解决这一问题有时还是相当困难的，
蒙特卡罗方法则不
存在这一问题。
蒙特卡罗方法的费用：根据上述误差公式，若问题所要求的误差为，置信水平为1    -
，X是按照正态
积分表由置信水平所确定的，则应有
   
    即子样容量N必须满足
   
    进一步假设每计算一次无偏统计量所需要的费用是C，则蒙特卡罗方法的总费用自然
应该是
   
    就是说，蒙特卡罗方法的费用与方法所确定的无偏统计量的方差及其费用C的乘积C成
正比。
    上述结果表明，在相同条件下，即在相同误差和置信水平要求下，一个蒙特卡罗方法
的好坏完全取决于
C的值的大小，其值越小相应的方法越好。
蒙特卡罗方法的收敛速度与问题的维数无关：由蒙特卡罗方法的误差公式看出，在置信水
平一定的情况
下，蒙特卡罗方法的误差除了与问题所确定的无偏统计量的方差有关以外，只取决于子样
的容量N，而与
无偏统计量g()中的是多少维空间的点无关。例如，如下的S重积分计算问题：
   
   蒙特卡罗方法的近似估计是
   
    其中，为S维方体内均匀抽样确定的点。由蒙特卡罗方法的误差公式不难看出，
如果无偏统计量g(x)的方差不变，则除了由于维数的变化可能引起每次观察费用作较小的
变化之外，蒙特卡
罗方法的误差是与维数S无关的，或者说，蒙特卡罗方法的收敛速度与问题的维数无关，
显然，这一特点是
其他计算方法所不常具有的。
蒙特卡罗方法受问题的条件影响不大：蒙特卡罗方法的另一个基本特点是，它受问题条件
限制的影响不
大。例如，对于上述积分，如果加上一个积分区域V的限制：
   
    则无论S维单位方体内的区域V如何特殊，我们都可以给出类似的近似估计如下：
   
    显然，这也是其他计算方法所不常具有的。 蒙特卡罗方法不必进行离散化处理：蒙
特卡罗方法的再一
个特点是，它不象其他计算方法那样对问题一定要进行离散化处理。明确些说，其他计算
方法需要事先确定
网格点，一旦网格点确定后，问题所关心的就是这些网格点上的值，而与其他点上的值无
关。蒙特卡罗方法
完全不需要这些。例如，在化工问题中的如下齐次积分方程本征值：
   
    其中和x为三维空间的点，V一般都比较复杂。对此问题，一般计算方法总是需要在积
分区域V上划分
网格点(这一步用计算机实现是比较困难的)，实现离散化，使问题所关心的只是这些网格
点上的值(需要占
用计算机大量存贮单元)。蒙特卡罗方法则不需要这样做。很明显，蒙特卡罗方法的这一
优点，不仅使蒙特
卡罗方法更加适于解决那些高维问题，精确度进一步增加，节省计算机大量存贮单元，而
且给建立蒙特卡
罗方法通用程序带来了很大方便。
蒙特卡罗方法具有直接解决问题的能力：蒙特卡罗方法还有一个特点是，对于那些本身具
有统计性质的
所谓非确定性问题，不需要象常规方法那样首先将它转化成为确定性问题，如转化成某方
程的解，然后再
通过解决确定性问题得到问题的答案。例如，粒予进入屏蔽层与介质发生多次碰撞后可能
被吸收，也可能
反射回来或穿过屏蔽层，问题所关心的常是穿透概率。对于这样一个非确定性问题，常规
方法是首先给出
它所满足的方程，实际上是一个积分微分方程，即转化成为确定性问题，然后用某种计算
方法解这一确定
性问题，进而得到问题的解答。蒙特卡罗方法不需要将非确定性问题转化成为确定性问题
，可以直接从非
确定性问题出发，通过模拟原问题的实际过程得到问题的解决。在由非确定性问题转化成
确定性问题和用
计算方法解这一确定性问题过程中，常常需要作很多近似，蒙特卡罗方法则不需要，因此
，有人称蒙特卡
罗方法为精确的方法，并常依此为标准衡量其他方法的近似是否是合理的。
蒙特卡罗方法具有同时计算多个方案的能力：蒙特卡罗方法的再一个特点是，对于那些需
要计算多个方
案的问题，有时不需要象常规方法那样逐个方案遂个计算，而可以同时计算所有的方案，
其全部计算量几
乎与计算一个方案的计算量相当。例如，上述穿透概率计算问题，当屏蔽层为均匀介质平
几何情况，要计
算I种厚度的穿透概率时，用蒙特卡罗方法只须计算最厚一种情况，其他厚度的穿透概率
在计算最厚一种情
况时稍加处理便可同时得到。详细些说，在计算最厚一种屏蔽层的穿透概率时，粒子在屏
蔽层内随机游动，
对于各种厚度的屏蔽层，粒子是离开，或离开再进入，也可能再离开等等情况是请清楚楚
的。设第一次离
开各种厚度的粒子数依次为，观察的粒子总数为N，则，便分别近似等于各种
厚度屏蔽层的穿透概率。
蒙特卡罗方法的缺点：蒙特卡罗方法也存在一些缺点，这些缺点主要有，对于维数少的问
题，一般是二维
和二维以下问题，它不如其他计算方法好；对于大的几何系统或小概率事件的计算问题，
它的计算结果有
时比真值偏低；误差是概率误差，而不是一般意义下的误差。
3. 伪随机数
    我们已经看到，由具有已知分布的总体中产生简单子样，在蒙特卡罗方法的问题中占
有非常重要的地
位。总体和子样的关系属于一般和若干个别的关系，或者说属于共性和若干个性的关系。
由具有已知分布
的总体中产生简单子样，就是由简单子样中的诸个性近似地反映总体的共性。   
    在连续型分布中，最简单而又最基本的一个分布是单位均匀分布。在蒙特卡罗方法中
常把该分布的简单
子样作为基本量来看待，而由其他分布中产生简单子样时把它看成是已知的。因此，它们
虽然都属于由具有
已知分布的总体中产生简单子样的问题，但就产生的方法而言，它们之间却有着本质上的
差别。
    均匀分布也常称为矩形分布，其中最基本的是单位均匀分布。单位均匀分布的密度函
数如下：
     f(x) = 1，0≤x≤1。
    由具有单位均匀分布的总体中产生的简单子样：，简称为随机数序列，其中的每一个
体称为
随机数。随机数在蒙特卡罗方法中占有极其重要的地位。根据随机数的定义，随机数序列
是相互
独立的具有相同单位均匀分布的随机变量序列。随机数具有非常重要的性质：对于任意自
然数S，由S个随机
数所组成的S维空间上的点在S维空间的单位方体上均匀分布，亦即对于任意的，
0≤≤1，i = 1，…，S，如下等式成立：
   
    反之，如果随机变量序列，对于任意自然数S，由S个元素所组成的   
维空间上的点
在上均匀分布，则随机变数序列是随机数序列。
    在机器中若没有随机数发生器，可以用下面的子程序RAN产生随机数(在主程序中要先
调用INITRAN置初
值)。
   
4. 由已知分布的随机抽样
    由已知分布的随机抽样指的是由已知分布的总体中产生简单子样。用F(x)表示已知分
布的分布函数，
，表示由总体F(x)中产生的容量为N的简单子样。按照简单子样的定义，
是用相互独立的抽样方法产生的，具有相同的分布F(x)。由已知分布的随机抽样简称为随
机抽样。用表
示由已知分布F(x)所产生的简单子样中的个体。对于连续型分布常用密度函数f(x)表示总

体的已知分布。用表示由已知分布f(x)所产生的简单子样中的个体。
    抽样方法有好几种，如直接抽样方法、舍选抽样方法、复合抽样方法、复合舍选抽样
方法等。然而，
由于在实际问题中随机变量所服从的分布是千差万别的，用这些方法实现随机抽样有时还
会遇到困难。概
括起来说，可能遇到的困难主要有以下两点：一点是有的分布用上述抽样方法虽然可以实
现，但抽样效率
很低；另一点是有的分布用上述方法虽然可以实现，抽样效率也不低，但运算量太大。由
于上述原因，常
采用某种近似抽样方法。



5. 系综的平均观测量

6. 线型高分子链构象模拟示例
    线型高分子是由数目很大的结构单元连接而成的长链分子，例如分子量为28,000的聚
乙烯(PE)，结构单
元数 = 1000。由于热运动及分子中的单键可以旋转而产生许多不同的构象。当分子量确
定后，随着构象的
变化，分子的尺寸也在变化。我们用分子末端距（线型高分子链两端的直线距离，h）的
均方根来作为描述
分子尺寸的参数，通常称为根均方末端距。一般情况下h应该是分子内化学键向量和的模
长（e是单位向量，
n和l分别为键数和键长）
   
    高分子链完全伸直时，主链为锯齿状，对聚乙烯有，l =    0.154nm。
    根据高分子构象的一些理论模型可以导出分子末端距。例如：

    柔性链模型的模拟则是在以链段长为半径的球面上取均匀分布的随机点，作为下一链
段的起点。同样
在n步以后，计算该链的末端距平方。统计了N条高分子链以后，计算出均方末端距。
7. 固体反应中成核过程模拟示例（化学物理学报，7/4 (1994)，118）
    固体反应与在流体中进行的反应有很大差别，一般分为成核与生长两个阶段。
   在较低温度下，分散在母相中的产物以替代缺陷形式占据晶格点，构成二元物系。在
成对最近邻相互作
用近似下，各向同性互作用系的哈密顿可以表示为
     
    式中A和B分别是缺陷在格点上的形成能和成对能；当i格点被产物缺陷占据时，占有
数= 1，否则，被
反应物占据时，占有数 = 0；表示对所有格点位置求和， 表示求和仅对i的最近邻进行。
因而
分别是缺陷总数和缺陷对数。哈密顿第二项前的负号表示产物缺陷成对使体系
能量降低，这是产物相能成核从反应物相中分离出来的原因。

wd840611

<>谢了</P>

hwlls

<>强</P>

wangguanjin

<>不错</P>

bob95

<>多谢了</P>

wujinglu

<>thank you very much!</P>

xiaoshiyilang

<>谢谢了!有兴趣的就看一看吧!参赛的就更应该学了</P>

datouyxa

公式无法显示，倘若直接上传文件就好了

天山雪狼

公式都没有！瑕疵呀！

天山雪狼

<>可不可以发个完整的出来！谢谢！</P><>或者我的邮箱：<a href="mailtnykie@163.com" target="_blank" >nykie@163.com</A>，万分感谢！</P>