《MonteCarlo方法》PPT课件.ppt

Monte Carlo 方法 黄世萍 起源 这一方法源于美国在第二次世界大战进研制原子弹的 曼哈顿计划 。 Monte Carlo方法创始人主要是这四位： Stanislaw Marcin Ulam, Enrico Fermi, John von Neumann和 Nicholas Metropolis。 Monte Carlo 方法基础 Monte Carlo（ MC） 方法是在简单的理论准则基础上 （ 如简单 的物质与物质以及物质与环境相互作用 ） ， 采用反复随机抽样 的手段 ， 解决复杂系统的问题 。 该方法采用随机抽样的手法 ， 可以模拟对象的概率与统计的问题 。 通过设计适当的概率模型 ， 该方法还可以解决确定性问题 ， 如定积分等 。 随着计算机的迅 速发展 ， MC 方法已在应用物理 、 原子能 、 固体物理 、 化学 、 材料 、 生物 、 生态学 、 社会学以及经济学等领域得到了广泛的 应用 。 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 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 geometrical analysis; kinetic Monte Carlo: simulate processes using scaling arguments to establish timescales or by introducing stochastic effects into molecular dynamics. Monte Carlo方法的基本思想 Monte Carlo 方法的基本思想是 : 为了求解某个问题 ， 建立一个恰 当的概率模型或随机过程 ， 使得其参量 （ 如事件的概率 、 随机变 量的数学期望等 ） 等于所求问题的解 ， 然后对模型或过程进行反 复多次的随机抽样试验 ， 并对结果进行统计分析 ， 最后计算所求 参量 ， 得到问题的近似解 。 Monte Carlo 方法是随机模拟方法 ； 它不仅限于模拟随机性问题 ， 还可 以解决确定性的数学问题 。 对随机性问题 ， 可以根据实际问题的概率法 则 ， 直接进行随机抽样试验 ， 即直接模拟方法 。 对于确定性问题采用 间接模拟方法 ， 即通过统计分析随机抽样的结果获得确定性问题的解 。 用 Monte Carlo 方法解决确定性的问题主要是在数学领域 ， 如计算重积 分 、 求逆矩阵 、 解线性代数方程组 、 解积分方程 、 解偏微分方程边界问题 和计算微分算子的特征值等 。 用 Monte Carlo 方法解决随机性问题则在 众多的科学及应用技术领域得到广泛的应用 ， 如中子在介质中的扩散问 题 、 库存问题 、 随机服务系统中的排队问题 、 动物的生态竞争 、 传染病的 蔓延等 。 简单的例子 对积分进行变换 ， 构造新的被积函数 g（ x） ， 使得该函数满足下列条件 ： g（ x）是连续随机变量 的概率密度函数。 定积分是概率积分 ， 其积分值等于概率 Pr （ a b） ， 即 这个步骤就是将一个 积分转化为一个概率模型的过程 ； 然后 ， 反复多次的 随机抽样试验 ，以抽样结果的统计平均作为索求概率的近似值 ， 从而求得 该积分。 具体试验步骤如下 ： （ 1） 产生服从给定分布函数 g（ x）的随机变量值 xi （ 2）检查 xi是否落入积分区域（ a x b） ， 如果满足条件 ， 则记录一次。 反复进行上述试验。 假设在 N 次试验后 ， xi落入积分区域的总次数为 m ， 那么 ， 积分值近似表示为 对于随机性问题 ， 可直接将实际的随机问题抽象为概率数学模型 ， 然后 与求解确定性问题一样进行抽样试验和统计计算。 Monte Carlo 方法解决实际问题的过程中 ， 主要有以下几个内容 建立简单而又便于实现的概率统计模型 ， 使所求的解正是该模型的某 一事件的概率或数学期望 ， 或该模型能够直接描述实际的随机过程 。 根据概率统计模型的特点和计算的需求 ， 改进模型 ， 以便减小方差和 减低费用 ， 提高计算效率 。 建立随机变量的抽样方法 ， 包括伪随机数和服从特定分布的随机变量 的产生方法 。 给出统计估计值及其方差或标准误差 。 One-Dimensional Integrals Methodical approaches rectangle rule, trapezoid rule, Simpsons rule Quadrature formula n uniformly separated points Sum areas of shapes approximating shape of curve () b a I f x dx Evaluating the general integral ()fx x 11 ( ) ( )nnii ii baI x f x f x n 假设对于下列积分 () b a I f x dx 首先，按照均匀分布，随即采样区域 a,b，记采样点为 Xi，每个采样点的 值 p(x)=1/(b-a)。利用 Monte Carlo算法，上述积分可以模拟为： 当 N趋近于无穷大的时候，我们认为两者的值是一致的，所以可以利用上 述离散的方式模拟这个积分 Monte Carlo Integration Stochastic approach Same quadrature formula, different selection of points 1 ()n i i baI f x n n points selected from uniform distribution p(x) x p(x) 我们在对于 a,b进行采样的时候 ， 完全没有必要进行均匀采样 ， 这样 做只是简单而已 ， 我们可以根据一种概率分布来进行采样 ， 从而使得 某些区域的采样密度更大 。 加入采样点的概率分布函数为 P(x)， 那么 我们可以利用如下公式计算定积分的值： 证明： 3 Monte Carlo方法的收敛性和基本特点 设所求的量 x是随机变量 的数学期望 E（ x） ， 那么 ， Monte Carlo 方法通常使用随机变量 的简单子样 ， ， ， N 的算术平均值 ， 即 作为所求量 x 的近似值。 由柯尔莫哥罗夫（ Kolmogorov）大数定理可知 ， 即当 N 充分大时 ， 有 成立的概率等于 ， 亦即可以用 N 作为所求量 x 的估计值。 根据中心极限定理 ， 如果随机变量 的标准差 不为零 ， 那么 Monte Carlo 方法的误差 为 式中 ， 为正态差 ， 是与置信水平有关的常量。 Monte Carlo 方法 的收敛速度的阶为 o（ N -1/2） ， 误差是由随机变量的标准差 s和抽样次 数 N 决定的。 提高精度一位数 ， 抽样次数要增加 倍 ； 减小随机变量的标准 差 ， 可以减小误差 。 Monte Carlo 方法具有以下四个重要特征 ： 由于 Monte Carlo 方法是通过大量简单的重复抽样来实现的 ， 因 此 ， 方法和程序的结构十分简单 。 收敛速度比较慢 ， 因此 ， 较适用于求解精度要求不高的问题 。 收敛速度与问题的维数无关 ， 因此 ， 较适用于求解多维问题 。 问题的求解过程取决于所构造的概率模型 ， 而受问题条件限制的 影响较小 ， 因此 ， 对各种问题的适应性很强 。 随机数的产生 1 随机数与伪随机数 Monte Carlo 方法的核心是 随机抽样 。 在该过程中往往需要各种各样分 布的随机变量其中最简单 、 最基本的是在 ， 区间上均匀分布的 随机变量 。 在该随机变量总体中抽取的子样 ， ， ， N 称为随 机数序列 ， 其中每个个体称为随机数 。 用数学的方法产生随机数是目前广泛使用的方法。 该方法的基本思想 是利用一种递推公式 ： 对于给定的初始值 ， 逐个地产生 ， ， 数学方法产生的随机数存在两个问题 ： 整个随机数序列是完全由递推函数形式和初始值唯一确定的 ， 严格 地说不满足随机数的相互独立的要求 。 存在周期现象 。 基于这两个原因 ， 将用数学方法所产生的随机数称为 伪随机数 。 伪 随机数的优点是适用于计算机 ， 产生速度快 ， 费用低廉 。 目前 ， 多数计算机均附带有 “ 随机数发生器 ” 。 选择递推函数必须注意以下几点 ： 随机性好 ； 在计算机上容易实现 ； 省时 ； 伪随机数的周期长。 2 伪随机数的产生方法 最基本的伪随机数是均匀分布的伪随机数。 该方法是首先给一个 r位的数 ， 取其中间的 r位数码作为第一个伪随机数 ， 然后将这个数平方 ， 构成一个新的 r位的数 ， 再取中间的 r位数作为第 二个伪随机数。 如此循环可得到一个伪随机数序列。 该方法的递推公式为 x表示对 x 取整 ， 运算 B（ Mod M）表示 B被 M 整除后的余数。 数列 i 是分布在 ，上的。 该方法由于 效率较低 ， 有时周期较短 ， 甚至会出现零 。 同余法 该方法也是由选定的初始值出发 ， 通过递推产生伪随机数序列。由于 该递推公式可写成数论中的同余式 ， 故称同余法。 该方法的递推公式 为 a ， c ， m分别称作倍数（ multiplier） 、增值（ Increment ）和模 （ modulus） ， 均为正整数 。 x 称为种子或初值 ， 也为正整数。 该方法所产生伪随机数的质量 ， 如周期的长度、独立性和均匀性 都与式中三个参数有关。 3 伪随机数的统计检验 伪随机数的特性好坏将直接影响到 Monte Carlo 的计算结果 ， 因此要对 所产生的伪随机数序列进行随机性检验 。 随机性检验主要包括 均匀性检 验 、 独立性检验 、 组合规律性检验和无连贯性检验 。 检验是伪随机数 检验最常用的方法 。 1. 均匀性 就是伪随机数列的 N 个数是否均匀分布在 ， 区间上。 若将 ，区间分成 k个相等的子区间（一般 k ， ， ） ， 若所得伪随机数在 ，区间 上是均匀分布的 ， 则虚假设 H 应为“每个伪随机数属于第 i组 的概率为 2. 频率检验 检验每组观测频数 ni与理论频数 mi N /k之间相差的显著性 3. 独立性 按先后顺序排列的 N 个伪随机数中 ， 每个数的出现是否与 其前后各个数独立无关。 对于两组伪随机数来说 ， 独立性 就是指它们不相关 。 4. 组合规律性 将 N 个伪随机数按一定的规律组合起来 ， 则各种组合的出现具 有一定的概率。 5. 无连贯性 将一次出现的 N 个伪随机数 ， 按其大小分为两类或 k类 ， 则各 类数的出现是否没有连贯现象。 确定性问题的 Monte Carlo 方法求解 Monte Carlo 方法所能解决的问题可以分为两大类 确定性问题 随机性问题 确定性问题主要包括： 求解线性和非线性方程组、逆矩阵、椭圆型差分方程的边值、积分方程 以及多重积分等。 用 Monte Carlo 方法求解确定性问题的基本思想是 ： 对于给定的确定性问 题 ， 设计一个 概率统计模型 （ 或随机过程模型 ） ， 然后采用一定的抽样 方法 ， 按照所设计的概率统计模型进行抽样 ， 最后把这个模型产生的一 个数字特征作为该确定性问题的近似解 。 蒲丰试验 年法国著名学者蒲丰（ Buffon）就发表了采用随机抽样法计算圆 周率的论文。 这就是蒲丰随机投针试验 ， 即著名的蒲丰问题。 蒲丰概率模型是 : 在平面上画相距均为 a的平行线束 ， 向该平面上随 机投置一枚长为 l 的针 。 为了避免针同时与两条平行线相交的复杂情 况 ， 设定 l a。 如图所示 ， M 为针的中点 ， y为针的中点 M 到与之最近的平行线的距离 ， 为针与平行线的夹角 。 显然 ， y a ， 。 该随机试验所有可能的集为 针与平行线相交的充分必要条件是 y lsin ， 针与平行线相交事件的集为 则针与平行线相交的概率为 由于 l和 a均为已知常数 ， 只要通过大量抽样试验求得该概率 p ， 由上式即 可算出圆周率。 设投针总次数为 N ， 其中针与平行线相交次数为 v ， 由贝 努里（ Bernoulli）定理可知 ， 当 N充分大时 ， 该事件出现的频率接近于其概 率 ， 即 这就是蒲丰试验求圆周率的过程。虽然 ， 该方法要想获得较高的精确度 ， 需要数以百万次的抽样试验 ， 效率很低 ， 但蒲丰试验具有 Monte Carlo 方 法解决确定性问题的基本思想。 A Classic Example The Calculation of 随机性问题的 Monte Carlo 模拟 该过程是先建立一个 随机过程模型 ， 使得该过程的随机变量的数学期望 等于所要求解确定问题的解。这种计算方法叫做间接模拟方法 。 另一方 面 ， 采用随机试验的方法直接模拟随机过程 ， 解决随机问题 ， 这就 是所谓直接 Monte Carlo 模拟。 如模拟布朗运动、扩散过程、有机高分 子形态、晶粒生长等随机问题的模拟。 1 随机行走（ random walk）模拟 随机行走是一种典型的简单抽样方法 ， 可用以模拟扩散 、溶液中长 而柔性的大分子的性质等 . 随机行走主要有三种类型 : 无限制的随机行走（ Unrestricted random walk ， RW） 不退行走（ Nonreversal random walk ， NRRW） 自回避行走（ Self-avoiding walk ， SAW） 无限制随机行走就是指 ， 某一个质点的每一次行走没有任何限制 ， 既与 前一次行走无关 ， 也与以前任何一步所到的位置无关 。 这种模型可以用于 模拟质点的扩散等过程 ， 但是 ， 不能用于模拟高分子的位形 。 因为 ， 用随机行走方法模拟高分子位形是用随机行走的轨迹代表高分子的位形 ， 行走过的位置代表的是构成分子的原子或官能团 ， 因此 ， 无限制随机行 走忽略了体斥效应 。 不退行走就是禁止在每一步行走后立即倒退 ， 可以解决刚走的一步与上一 步重叠的问题 。 但不退行走没有完全解决高分子的体斥效应问题 。 自回避行走就是所有已走过的位置不能再走 ， 这样就完全解决了体斥效应 问题 。 二维方格子上的三种随机行走的示意图。 可以用四个矢量记述从某个节点 向邻近节点的行走方向（设格子间距为） ： N 步无限制随机行走的算法如下 ： 取 r （ ， ） （ 坐标原点 ） ， 并令 k ； 取一个在 和 之间的随机整数 vk ； k k ， rk rk （ vk ） ； 若 k N ， 行走结束 ， 否则回到第 步 。 对于不退行走 ， 可选择的行走方向不再是 ， 而是 。 禁止在每一步 行走后立即倒退 ， 即第 k步的方向矢量不能与第 k 步的方向矢量相 逆。 由式（ ）可以看出方向矢量 （） ， （）互为逆方 向 ， （） ， （）互为逆方向。 2 Markov链 对于简单抽样 ， 每一次的抽样都是独立的 。 如上述的随机行走过 程 ， 每行走一步都与前一步无关 ， 更与初始位置无关 。 Sampling 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 the 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 ! () N N N N U rZNU d r U r e 0Ue Importance Sampling Put more quadrature points in regions where integral receives its greatest contributions Return to 1-dimensional example Most contribution from region near x = 1 Choose quadrature points not uniformly, but according to distribution (x) linear form is one possibility How to revise the integral to remove the bias? 1 2 0 3I x dx 2( ) 3f x x ( ) 2xx The 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 sampling 1 ()n ii i I f x x 1x 2x 3x nx 1 ()i i bax nx Greater more points smaller spacing 1 () () () n i ii x b a f xI nx choose x points according to 在正则系综中 ， 任意观察量 A（ x）的热平均为 重要性抽样： 在随机过程中 ， 选取一个随机抽样的分布 ， 使生 成的随机数满足选取分布形式 。 根据一定的分布 形式进行的随机抽样称为重要性抽样 。 重要抽样 Monte Carlo 方法的实质是每次抽样试验不是完全独 立的 ， 而是与前一次或者与以前的所有抽样结果具有一定的 概率关系 ， 如不退随机行走和自回避随机行走 。 设一个系统的状态序列 （ 随机变量序列 ） 为 x ， x ， ， xn ， ， 如果对于任何一个状态 xn只与前一个状态 xn 有 关 ， 而与初始状态无关 ， 即状态 n的概率为 则称此序列为 Markov 链。 Markov 链是一种随机行走状态 ， 从状态 i单步行走到状态 j 的概 率叫做转移概率 ， 或跃迁概率 ， 即 设所有可能的状态数为 N ， 由 pij 构成的 N N 矩阵叫做转移矩阵 p ， 该矩阵的每一行元素的和等于 。 Markov 链的重要性质是 ， 无论初始 状态如何 ， 最终状态（足够多的时间步长次数）会遵从某一个唯一的分 布 ， 该分布叫做极限分布 xlim ， 即 也就是说 ， 极限状态乘转移概率后状态不再发生变化 ， 即系统达到一 个平衡状态。 因此 ， Markov 链在平衡态 Monte Carlo 模拟中具有重要 的意义。 3 Metropolis Monte Carlo法 Monte Carlo 方法主要分为简单随机抽样方法和重要随机抽样方法 。 简单抽样就是以平均分布进行抽样 ， 每次抽样是完全独立的 。 正 如前面关于积分问题中所述 ， 很多问题难以用简单抽样方法解决 ， 而重要随机抽样能够获得很好结果 。 Metropolis Monte Carlo 方法是一种重要随机抽样方法。 Metropolis 等人 提出了一种基于建立一个 Markov 过程的方法 。 该方法的实质是 ， 系统 的各状态不是彼此独立无关地选取 ，而是建造一个 Markov 过程 ， 过程 中每一个状态 xi 是由前一个状态 xi 通过一个适当的跃迁概率 W （ xi x i ）得到的 ， 并且 ， 该概率能够使得在 M 的极限下 ， Markov 过程产生的状态的分布函数 P（ xi ）趋于所要的平衡分布 ， 即 满足上述要求的充分条件为 也就是说 ， 两个状态正向与反向的跃迁概率之比只依赖于两者的 能量差 dH H （ xj ） H（ xi ） 。 但是满足该条件的跃迁概率 W 的形式并不是唯一 的 。 通常采用以下两种形式 ： ts 可取为 Metropolis Monte Carlo 方法的具体步骤如下 : Metropolis Monte Carlo （） 建立体系状态与能量的关系模型。 （） 由初始状态出发 ， 通过简单抽样设立新状态。 （） 根据新旧状态的哈密顿量 dH ， 判断新状态的舍选 ， 判断 舍选有以下 种情况 ： dH ， 接受新状态 ， 并在该状态基础上 ， 进一步进 行步骤（） ； dH ， 不是直接否决 ， 而是进一步判断 ， 抽取一个 随机数 如果新状态被拒绝 ， 则把原来的状态作为新状态 ， 重复进行步骤 （） ， 并记录一次。 如果系统的粒子数为 M ， 每次新状态的抽样均随机抽选一个 粒子 ， 并不是每个粒子逐一地进行 。 只要伪随机数的质量足 够高 ， 各粒子被抽样的概率是均等的 。 当抽样次数达到系统 粒子总数 M 时 该过程叫做一个 Monte Carlo step MCS Metropolis 方法在状态抽样时虽然采用的是简单抽样 ， 但是通过新旧状 态的能量判断 ， 实现新状态的舍选 ， 建立 Markov 过程 。 对于恒定组 成的正则和微正则系综 ， 系统的能量用哈密顿量表示 。 对于变化组成 的巨正则系综 ， 随机选取一个粒子并通过改变粒子的种类得到一个新的 组态 ， 系统的能量用混合能及混合物化学位之和表示 。