Monte Carlo 方法PPT课件

1、起源 这一方法源于美国在第二次世界大战进研制原子弹的曼哈顿计划。Monte Carlo方法创始人主要是这四位：Stanislaw Marcin Ulam, Enrico Fermi, John von Neumann和 Nicholas Metropolis。 第1页/共48页 Monte Carlo 方法基础Monte Carlo（MC）方法是在简单的理论准则基础上（如简单的物质与物质以及物质与环境相互作用） ， 采用反复随机抽样的手段 ， 解决复杂系统的问题 。 该方法采用随机抽样的手法 ，可以模拟对象的概率与统计的问题。通过设计适当的概率模型 ， 该方法还可以解决确定性问题 ， 如定积分

2、等。 随着计算机的迅速发展 ， MC 方法已在应用物理、原子能、固体物理 、化学、材料、生物、生态学、社会学以及经济学等领域得到了广泛的应用。第2页/共48页MC meothds classical Monte Carlo: samples are drawn from a probability distribution, often the classical Boltzmann distribution, to obtain thermodynamic properties or minimum-energy structures;quantum Monte Carlo: random

3、walks are used to compute quantum-mechanical energies and wavefunctions, often to solve electronic structure problems, using Schrdingers equation as a formal starting point; volumetric Monte Carlo: random number generators are used to generate volumes per atom or to perform other types of geometrica

4、l analysis;kinetic Monte Carlo: simulate processes using scaling arguments to establish timescales or by introducing stochastic effects into molecular dynamics.第3页/共48页Monte Carlo方法的基本思想Monte Carlo 方法的基本思想是: 为了求解某个问题 ， 建立一个恰当的概率模型或随机过程 ， 使得其参量（如事件的概率、随机变量的数学期望等）等于所求问题的解 ， 然后对模型或过程进行反复多次的随机抽样试验 ， 并对结

5、果进行统计分析 ， 最后计算所求参量 ， 得到问题的近似解。第4页/共48页Monte Carlo 方法是随机模拟方法 ； 它不仅限于模拟随机性问题 ， 还可以解决确定性的数学问题。 对随机性问题 ， 可以根据实际问题的概率法则 ， 直接进行随机抽样试验 ， 即直接模拟方法。 对于确定性问题采用间接模拟方法 ， 即通过统计分析随机抽样的结果获得确定性问题的解。用 Monte Carlo 方法解决确定性的问题主要是在数学领域 ， 如计算重积分、求逆矩阵、解线性代数方程组、解积分方程、解偏微分方程边界问题和计算微分算子的特征值等。 用 Monte Carlo 方法解决随机性问题则在众多的科学及应用

6、技术领域得到广泛的应用 ， 如中子在介质中的扩散问题、库存问题、随机服务系统中的排队问题、动物的生态竞争、传染病的蔓延等。第5页/共48页简单的例子对积分进行变换 ， 构造新的被积函数 g（x） ， 使得该函数满足下列条件 ：g（x）是连续随机变量 的概率密度函数。 定积分是概率积分 ， 其积分值等于概率 Pr （a b） ， 即这个步骤就是将一个积分转化为一个概率模型的过程 ； 然后 ， 反复多次的随机抽样试验 ，以抽样结果的统计平均作为索求概率的近似值 ， 从而求得该积分。第6页/共48页具体试验步骤如下 ：（1） 产生服从给定分布函数 g（x）的随机变量值 xi（2）检查 xi是否落入积

7、分区域（a x b） ， 如果满足条件 ， 则记录一次。反复进行上述试验。 假设在 N 次试验后 ， xi落入积分区域的总次数为 m ， 那么 ， 积分值近似表示为对于随机性问题 ， 可直接将实际的随机问题抽象为概率数学模型 ， 然后与求解确定性问题一样进行抽样试验和统计计算。第7页/共48页Monte Carlo 方法解决实际问题的过程中 ， 主要有以下几个内容 建立简单而又便于实现的概率统计模型 ， 使所求的解正是该模型的某 一事件的概率或数学期望 ， 或该模型能够直接描述实际的随机过程。 根据概率统计模型的特点和计算的需求 ， 改进模型 ， 以便减小方差和减低费用 ， 提高计算效率。 建

8、立随机变量的抽样方法 ， 包括伪随机数和服从特定分布的随机变量的产生方法。 给出统计估计值及其方差或标准误差。第8页/共48页One-Dimensional Integrals Methodical approaches rectangle rule, trapezoid rule, Simpsons rule Quadrature formulan uniformly separated pointsSum areas of shapes approximating shape of curve( )baIf x dxEvaluating the general integral( )f x

9、x11()()nniiiibaIxf xf xn 第9页/共48页假设对于下列积分 ( )baIf x dx首先，按照均匀分布，随即采样区域a,b，记采样点为Xi，每个采样点的值p(x)=1/(b-a)。利用Monte Carlo算法，上述积分可以模拟为： 当N趋近于无穷大的时候，我们认为两者的值是一致的，所以可以利用上述离散的方式模拟这个积分 第10页/共48页Monte Carlo Integration Stochastic approach Same quadrature formula, different selection of points1()niibaIf xnn poin

10、ts selected from uniform distribution p(x)xp(x)第11页/共48页我们在对于a,b进行采样的时候，完全没有必要进行均匀采样，这样做只是简单而已，我们可以根据一种概率分布来进行采样，从而使得某些区域的采样密度更大。加入采样点的概率分布函数为P(x)，那么我们可以利用如下公式计算定积分的值： 证明：第12页/共48页3 Monte Carlo方法的收敛性和基本特点设所求的量 x是随机变量的数学期望 E（x） ， 那么 ， Monte Carlo 方法通常使用随机变量的简单子样， ， ，N 的算术平均值 ， 即作为所求量 x 的近似值。 由柯尔莫哥罗夫（

11、Kolmogorov）大数定理可知 ，即当 N 充分大时 ， 有成立的概率等于 ， 亦即可以用 N 作为所求量 x 的估计值。第13页/共48页根据中心极限定理 ， 如果随机变量 的标准差 不为零 ， 那么 Monte Carlo 方法的误差为式中 ， 为正态差 ， 是与置信水平有关的常量。 Monte Carlo 方法的收敛速度的阶为 o（N -1/2） ， 误差是由随机变量的标准差 s和抽样次数 N 决定的。提高精度一位数 ， 抽样次数要增加 倍 ； 减小随机变量的标准差 ， 可以减小误差 。第14页/共48页Monte Carlo 方法具有以下四个重要特征 ： 由于 Monte Carl

12、o 方法是通过大量简单的重复抽样来实现的 ， 因此 ， 方法和程序的结构十分简单 。 收敛速度比较慢 ， 因此 ， 较适用于求解精度要求不高的问题。 收敛速度与问题的维数无关 ， 因此 ， 较适用于求解多维问题。 问题的求解过程取决于所构造的概率模型 ， 而受问题条件限制的影响较小 ， 因此 ， 对各种问题的适应性很强。第15页/共48页随机数的产生1 随机数与伪随机数Monte Carlo 方法的核心是随机抽样。 在该过程中往往需要各种各样分布的随机变量其中最简单、最基本的是在 ，区间上均匀分布的随机变量。 在该随机变量总体中抽取的子样 ， ， ，N 称为随机数序列 ， 其中每个个体称为随机

13、数。用数学的方法产生随机数是目前广泛使用的方法。 该方法的基本思想是利用一种递推公式 ：对于给定的初始值 ， 逐个地产生 ， ， 第16页/共48页数学方法产生的随机数存在两个问题 ： 整个随机数序列是完全由递推函数形式和初始值唯一确定的 ， 严格地说不满足随机数的相互独立的要求。 存在周期现象。基于这两个原因 ， 将用数学方法所产生的随机数称为伪随机数。 伪随机数的优点是适用于计算机 ， 产生速度快 ， 费用低廉。 目前 ， 多数计算机均附带有“随机数发生器” 。第17页/共48页选择递推函数必须注意以下几点 ： 随机性好 ； 在计算机上容易实现 ； 省时 ； 伪随机数的周期长。第18页/共

14、48页2 伪随机数的产生方法最基本的伪随机数是均匀分布的伪随机数。该方法是首先给一个 r位的数 ， 取其中间的 r位数码作为第一个伪随机数 ， 然后将这个数平方 ， 构成一个新的 r位的数 ， 再取中间的 r位数作为第二个伪随机数。 如此循环可得到一个伪随机数序列。 该方法的递推公式为x表示对 x 取整 ， 运算 B（Mod M）表示 B被 M 整除后的余数。 数列i 是分布在 ，上的。 该方法由于效率较低 ， 有时周期较短 ， 甚至会出现零 。第19页/共48页同余法该方法也是由选定的初始值出发 ， 通过递推产生伪随机数序列。由于该递推公式可写成数论中的同余式 ， 故称同余法。 该方法的递推

15、公式为a ， c ， m分别称作倍数（multiplier） 、增值（Increment ）和模（modulus） ， 均为正整数 。 x 称为种子或初值 ， 也为正整数。 该方法所产生伪随机数的质量 ， 如周期的长度、独立性和均匀性都与式中三个参数有关。第20页/共48页第21页/共48页3 伪随机数的统计检验伪随机数的特性好坏将直接影响到 Monte Carlo 的计算结果 ， 因此要对所产生的伪随机数序列进行随机性检验。 随机性检验主要包括均匀性检验 、独立性检验、组合规律性检验和无连贯性检验。 检验是伪随机数检验最常用的方法。1. 均匀性就是伪随机数列的 N 个数是否均匀分布在 ，区间

16、上。 若将 ，区间分成 k个相等的子区间（一般 k ， ， ） ， 若所得伪随机数在 ，区间上是均匀分布的 ， 则虚假设H 应为“每个伪随机数属于第 i组的概率为第22页/共48页2. 频率检验 检验每组观测频数 ni与理论频数mi N /k之间相差的显著性3. 独立性按先后顺序排列的 N 个伪随机数中 ， 每个数的出现是否与其前后各个数独立无关。 对于两组伪随机数来说 ， 独立性就是指它们不相关 。4. 组合规律性将 N 个伪随机数按一定的规律组合起来 ， 则各种组合的出现具有一定的概率。5. 无连贯性将一次出现的 N 个伪随机数 ， 按其大小分为两类或 k类 ， 则各类数的出现是否没有连贯

17、现象。第23页/共48页确定性问题的 Monte Carlo 方法求解Monte Carlo 方法所能解决的问题可以分为两大类确定性问题随机性问题确定性问题主要包括： 求解线性和非线性方程组、逆矩阵、椭圆型差分方程的边值、积分方程以及多重积分等。用 Monte Carlo 方法求解确定性问题的基本思想是 ： 对于给定的确定性问题 ， 设计一个概率统计模型（或随机过程模型） ， 然后采用一定的抽样方法 ， 按照所设计的概率统计模型进行抽样 ， 最后把这个模型产生的一个数字特征作为该确定性问题的近似解。第24页/共48页蒲丰试验 年法国著名学者蒲丰（Buffon）就发表了采用随机抽样法计算圆周率的

18、论文。 这就是蒲丰随机投针试验 ， 即著名的蒲丰问题。蒲丰概率模型是 : 在平面上画相距均为a的平行线束 ， 向该平面上随机投置一枚长为 l 的针。 为了避免针同时与两条平行线相交的复杂情况 ， 设定 l a。 如图所示 ， M 为针的中点 ， y为针的中点M 到与之最近的平行线的距离 ， 为针与平行线的夹角。 显然 ， y a ， 。 该随机试验所有可能的集为第25页/共48页针与平行线相交的充分必要条件是 y lsin ，针与平行线相交事件的集为则针与平行线相交的概率为由于 l和 a均为已知常数 ， 只要通过大量抽样试验求得该概率 p ， 由上式即可算出圆周率。 设投针总次数为 N ， 其

19、中针与平行线相交次数为 v ， 由贝努里（Bernoulli）定理可知 ， 当 N充分大时 ， 该事件出现的频率接近于其概率 ， 即这就是蒲丰试验求圆周率的过程。虽然 ， 该方法要想获得较高的精确度 ， 需要数以百万次的抽样试验 ， 效率很低 ， 但蒲丰试验具有 Monte Carlo 方法解决确定性问题的基本思想。第26页/共48页A Classic Example The Calculation of 第27页/共48页第28页/共48页随机性问题的 Monte Carlo 模拟该过程是先建立一个随机过程模型 ， 使得该过程的随机变量的数学期望等于所要求解确定问题的解。这种计算方法叫做间接

20、模拟方法。 另一方面 ， 采用随机试验的方法直接模拟随机过程 ， 解决随机问题 ， 这就是所谓直接 Monte Carlo 模拟。 如模拟布朗运动、扩散过程、有机高分子形态、晶粒生长等随机问题的模拟。第29页/共48页1 随机行走（random walk）模拟随机行走是一种典型的简单抽样方法 ， 可用以模拟扩散 、溶液中长而柔性的大分子的性质等.随机行走主要有三种类型:无限制的随机行走（Unrestricted random walk ， RW）不退行走（Nonreversal random walk ， NRRW）自回避行走（Self-avoiding walk ， SAW）第30页/共48

21、页无限制随机行走就是指 ， 某一个质点的每一次行走没有任何限制 ， 既与前一次行走无关 ， 也与以前任何一步所到的位置无关。 这种模型可以用于模拟质点的扩散等过程 ， 但是 ， 不能用于模拟高分子的位形。 因为 ， 用随机行走方法模拟高分子位形是用随机行走的轨迹代表高分子的位形 ， 行走过的位置代表的是构成分子的原子或官能团 ， 因此 ， 无限制随机行走忽略了体斥效应。不退行走就是禁止在每一步行走后立即倒退 ， 可以解决刚走的一步与上一步重叠的问题。 但不退行走没有完全解决高分子的体斥效应问题。自回避行走就是所有已走过的位置不能再走 ， 这样就完全解决了体斥效应问题。第31页/共48页二维方格

22、子上的三种随机行走的示意图。 可以用四个矢量记述从某个节点向邻近节点的行走方向（设格子间距为） ：第32页/共48页N 步无限制随机行走的算法如下 ： 取 r （ ， ）（坐标原点） ， 并令 k ； 取一个在 和 之间的随机整数 vk ； k k ， rk rk （vk ） ； 若 k N ， 行走结束 ， 否则回到第 步。对于不退行走 ， 可选择的行走方向不再是 ， 而是 。 禁止在每一步行走后立即倒退 ， 即第k步的方向矢量不能与第 k 步的方向矢量相逆。 由式（ ）可以看出方向矢量 （） ， （）互为逆方向 ， （） ， （）互为逆方向。第33页/共48页2 Markov链对于简单抽样

23、 ， 每一次的抽样都是独立的。如上述的随机行走过程 ， 每行走一步都与前一步无关 ， 更与初始位置无关。Sampling第34页/共48页Integrate Over a Simple Shape? Statistical-mechanics integrals typically have significant contributions from miniscule regions of the integration space contributions come only when no spheres overlap e.g., 100 spheres at freezing t

24、he fraction is 10-260 Evaluation of integral is possible only if restricted to region important to integral must contend with complex shape of region MC methods highly suited to “importance sampling”()11!()NNNNU rZNUdr U re0Ue第35页/共48页Importance Sampling Put more quadrature points in regions where i

25、ntegral receives its greatest contributions Return to 1-dimensional example Most contribution from region near x = 1 Choose quadrature pointsnot uniformly, but accordingto distribution (x) linear form is one possibility How to revise the integral to remove the bias?1203Ix dx2( )3f xx( )2xx第36页/共48页T

26、he Importance-Sampled Integral Consider a rectangle-rule quadrature with unevenly spaced abscissas Spacing between points reciprocal of local number of points per unit length Importance-sampled rectangle rule Same formula for MC sampling1()niiiIf xx1x2x3xnx1()iibaxnxGreater more points smaller spaci

27、ng1( )()()niiixbaf xInxchoose x points according to 第37页/共48页在正则系综中 ， 任意观察量 A（x）的热平均为第38页/共48页重要性抽样：在随机过程中，选取一个随机抽样的分布，使生成的随机数满足选取分布形式。根据一定的分布形式进行的随机抽样称为重要性抽样。第39页/共48页重要抽样 Monte Carlo 方法的实质是每次抽样试验不是完全独立的 ， 而是与前一次或者与以前的所有抽样结果具有一定的概率关系 ， 如不退随机行走和自回避随机行走。设一个系统的状态序列（随机变量序列）为 x ，x ， ，xn ， ， 如果对于任何一个状态 x

28、n只与前一个状态 xn 有关 ， 而与初始状态无关 ， 即状态 n的概率为则称此序列为 Markov 链。第40页/共48页Markov 链是一种随机行走状态 ， 从状态 i单步行走到状态 j 的概率叫做转移概率 ， 或跃迁概率 ， 即设所有可能的状态数为 N ， 由 pij 构成的 N N 矩阵叫做转移矩阵 p ， 该矩阵的每一行元素的和等于 。 Markov 链的重要性质是 ， 无论初始状态如何 ， 最终状态（足够多的时间步长次数）会遵从某一个唯一的分布 ， 该分布叫做极限分布 xlim ， 即也就是说 ， 极限状态乘转移概率后状态不再发生变化 ， 即系统达到一个平衡状态。 因此 ，Markov 链在平衡态 Monte Carlo 模拟中具有重要的意义。第41页/共48页3 Metropolis Monte Carlo法Monte Carlo 方法主要分为简单随机抽样方法和重要随机抽样方法。 简单抽样就是以平均分布进行抽样 ， 每次抽样是完全独立的。正如前面关于积分问题中所述 ， 很多问题难以用简单抽样方法解决 ， 而重要随机抽样能够获得很好结果。Metropolis Mon