随机模拟与蒙特卡罗法.doc

随机模拟与蒙特卡罗法分子聚集体或宏观物系的结构与物理化学性质可以通过求相应物理量的系综平均值得到。但是，系统的状态数常常是一个近乎天文的数字。比如，对一个大小的一个二维点阵，即使每个点只能取两个状态，系统可能取的状态总数为。这样，求一个平均值，就是用速度为每秒数万亿次的巨型计算机，它也要算上至少年。显然，这是根本无法做到的。计算机随机模拟方法是用一小部分“代表”状态上的算术平均值近似系综平均值。 计算机随机模拟是在电子计算机上对随机现象进行模拟并进而得到问题解答的方法的简称。其中所说的随机现象可以是原问题本身所固有的，也可以是人为建立起来的与原问题答案有一定关系的某随机现象。计算机随机模拟方法亦称蒙特卡罗方法(Monte Carlo Method)，后一名称已被国内外广大学者所普遍采用。蒙特卡罗方法的程序大都为专门应用设计，一般的应用程序不会太复杂，可以自己编写。一个模拟退火算法程序ASA(Adaptive Simulated Annealing)可以在网上免费得到(，)，用于许多优化的目的。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，0x1。 由具有单位均匀分布的总体中产生的简单子样：，简称为随机数序列，其中的每一个体称为随机数。随机数在蒙特卡罗方法中占有极其重要的地位。根据随机数的定义，随机数序列是相互独立的具有相同单位均匀分布的随机变量序列。随机数具有非常重要的性质：对于任意自然数S，由S个随机数所组成的S维空间上的点在S维空间的单位方体上均匀分布，亦即对于任意的，01，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. 固体