蒙特卡罗方法课件(PPT 56页)

1、1目 录 第一章 蒙特卡罗方法概述第二章 随机数的产生第三章 EM算法和MCMC方法参考书 : 茆诗松等, 高等数理统计(第6章), 高等教育出版社,1998;2.徐钟济，蒙特卡罗方法，上海科学技术出版社第1页，共56页。2第一章 蒙特卡罗方法概述 蒙特卡罗方法又称随机抽样技巧或统计试验方法。 蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它以概率统计理论为基础。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。第2页，共56页。31.蒙特卡罗方法的基本思想 理论基础：大数定律；中心极限定理； F(X)U(0

2、,1)。基本思想：1.当所求问题的解是某个事件的概率，或者是某个随机变量的期望，或与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或该随机变量若干个观察值的算术平均值，根据大数定律得到问题的解；2. 要生成分布函数为F(x)的随机数，可先生成U(0,1)随机数F，则可得到随机数X=F-1(F) 。第3页，共56页。4例（利用MC进行欧式期权定价）设股票价格St服从风险中性测度下的几何Brown运动：其离散化形式为根据金融工程理论，设现在股票价格为S0，T时刻到期（单位天），敲定价为K的欧式看涨期权的价格为MC方案：按照（1）递推产生n条风险中性测度下的轨道，提取出ST (

3、n)；（2）第4页，共56页。52. 蒙特卡罗方法的误差 根据中心极限定理如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2 ，即则当N充分大时，有如下的近似式它表明，误差收敛速度的阶为 以概率1-成立。第5页，共56页。6通常，蒙特卡罗方法的误差定义为关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出 。 第6页，共56页。7减小方差的各种技巧 显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精

4、度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。降低方差的各种技巧，引起了人们的普遍注意。一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为 其中c是观察一个子样的平均费用。第7页，共56页。8蒙特卡罗方法的特点优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。受几何条件限制小。收敛速度与问题的维数无关。误差容易确定。程序结构简单，易于实现。 缺点收敛速度慢。误差具有概率性。第8页，共

5、56页。9第二章 随机数的产生2.1 逆变换法设随机变量X的分布函数为F(x),定义定理2.1 设随机变量U服从U(0,1)分布，则的分布函数为F(x).由定理2.1,要生成分布函数为F(x)的随机数，可先生成U(0,1)随机数U，则可得到随机数X=F-1(U) 第9页，共56页。102.2 合成法如果X的密度函数p(x)难于抽样,而X关于Y的条件密度函数p(x|y)以及Y的密度函数g(y)均易于抽样，则X 的随机数可如下产生：Step1 由Y的分布g(y)抽取y;Step2 由X关于Y的条件密度函数p(x|y)抽取x.例2.1 设X的密度函数为由合成法，X的随机数可如下抽取：1）取uU(0,

6、1); 2)取 ,确定i，使3) 由pi(x)抽取x.第10页，共56页。112.3 筛选抽样 当p(x)难以直接抽样时，如果可以将p(x) 表示成p(x)=ch(x)g(x),其中h(.)是一密度函数且易于抽样，而0g(y),回到1）上述方法就是筛选抽样法，它是一种非常重要的抽样方法，可解决许多难以直接抽样的分布的抽样问题。第11页，共56页。12h(x)的的选取有多种方法。一种直观的方法是：如果存在一个函数M(x)，满足p(x)M(x),且令h(x)=M(x)/c, 若h(x)易于抽样，则筛选抽样变为1）由U(0,1)抽取u，由h(y)抽取y；2）如果up(y)/M(y),则x=y停止；3

7、）如果u p(y)/M(y),回到1）。筛选抽样的理论依据如下：定理 设X的密度函数为p(x),且p(x)=ch(x)g(x)，其中01时，如果 ，则x=y， 否则转到1）；第14页，共56页。152.4 随机向量的抽样法设X1,Xk的联合概率密度为定理2.4 设U1,Uk是独立同分布的U(0,1)变量， X1,Xk是方程的解，其中 是对应于 的分布函数，则X1,Xk的分布为(2.4).(2.4)(2.5)第15页，共56页。16随机向量的逆变换抽样法：由U(0,1)分布独立地抽取u1,uk;用方程（2.5)解x1,xk例2.3 设X1,X2的联合密度函数为试生成X1,X2的随机数。解：第16

8、页，共56页。17相应的边际分布函数和条件分布函数分别为方程(2.5)变为此方程不易解，不妨交换两自变量的次序第17页，共56页。18相应的边际分布函数和条件分布函数分别为方程(2.5)变为对服从特定分布的随机向量有一些特殊的抽样方法。第18页，共56页。19例2.6 试生成k维正态分布 的随机数。解：注意到若 ，则存在下三角阵使其中C可由迭代实现：首先，由 ，有从而。因于是得依此类推，第19页，共56页。20一般迭代公式为至此，我们可以给出k维正态分布的抽样步骤：1)迭代计算 ；2)由N(0,1)分布独立抽取k个随机数 ；3)计算第20页，共56页。212.5 随机模拟计算2.5.1 随机投

9、点法考虑积分 ，设a,b有限，0f(x)M,令=(x,y):axb,0yM,并设(X,Y)是在上均匀分布的二维随机向量，其联合密度函数为则易见， 是中曲线f(x)下方面积。 假设我们向中投点，若点落在y=f(x)下方称为中的，则点中的概率为第21页，共56页。22若我们进行了n次投点，其中n0次中的，则可以得到一个估计不难看出， 是的无偏估计，且其方差为(2.5.1)第22页，共56页。232.5.2 样本均值法于是，积分注意到，若XU(a,b),则由大数定律，若 ，则MC方法为：1) 独立产生n个U(a,b)随机数2)按(2.5.2)估计。(2.5.2)第23页，共56页。24可证，在0f(

10、x)M条件下，2.5.3 降低方差的技术Monte Carlo 方法中一类重要的研究课题是考虑一些降低估计方差的技术。常用的方法有：重要抽样法，分层抽样法，关联抽样法等。一 重要抽样法由上节，样本平均法比投点法有效，将样本平均法做更一般的推广，设g(x)是(a,b)上的密度函数，改写第24页，共56页。25由大数定律，若 ，则MC方法为：1)选择适当的g(x),独立产生n个g(x)随机数2)由(2.5.3)估计。显然(2.5.3)第25页，共56页。26从理论上看，因，若f(x)0,取则有因为未知，这是作不到的，但它提示我们取g(x)与f(x)形状接近，应能降低方差。这就是重要抽样法的基本思想

11、。其方差与g(x)有关。问题变为，如何选择g(.)使估计的方差最小。第26页，共56页。27例2.5.1 分别用投点法，均值法，重要抽样法，求积分 ,比较各种方法的有效性。解 i)投点法1)产生随机数 2) 对每对 ，记 的次数为n0.则Gii)均值法1)产生随机数 2)第27页，共56页。28iii)重要抽样法由重要抽样法的思想，需选择一个与 相似的密度函数。由Taylor展开式 取1)产生随机数 2) 取则（数值计算）真值投点法均值法重要抽样法1.718281.87561.81451.7219模拟结果第28页，共56页。29二、分层抽样法另一种利用贡献率大小来降低估计方差的方法是分层抽样法

12、。它首先把样本空间D分成一些不交的小区间 ，然后在各小区间内的抽样数由其贡献大小决定。即，定义 ，则Di内的抽样数ni应与pi成正比。考虑积分将0,1分成m个小区间：则记 为第i个小区间的长度，i=1,m.在每个小区间上的积分值可用均值法估计出来，然后将其相加即可给出的一个估计。具体步骤为：第29页，共56页。301) 独立产生U(0,1)随机数2)计算3) 计算于是可得的估计为(2.5.4)易见， 是的无偏估计，其方差为(2.5.5)(2.5.6)第30页，共56页。31续例2.5.1 考察分层抽样法求积分 的方差。解：先将区间0,1划分成两个小区间0,0.5,0.5,1,则设一共抽n个随机

13、数，其中在0,0.5)上抽n1个，则使用分层抽样法求得 的方差为第31页，共56页。32对n1求导易知，在n固定下，当 时的方差最小，为如果我们将区间进行10等份，并确定出最优的抽样次数分配： ，则可得到分层抽样法估计的方差为.一般地，若诸 已知，在n固定下，当时，估计的方差最小，为第32页，共56页。33分层抽样法在实施上有两个主要问题，其一是怎样划分区间，简单而常用的方法是将区间等分；另一个问题是在区间划分好后如何确定抽样次数的分配。由于在实际中 总是未知的，因而前面最优分配的结论无法应用。即使如此，分层抽样法还是有其作用的。可以证明，即使取简单的分配也有事实上，取 ，代入(2.5.5)得

14、由Cauchy-Schwarz不等式，有据此，在(2.5.6)式两端各乘以 并相加得 于是第33页，共56页。34三、关联抽样法考虑积分差若用 估计，则其方差为显然，在 确定后， 正相关度越高，则 的方差越小。这便是关联抽样法的基本出发点。考虑用重要抽样法来估计I1,I2,即改写为产生n个U(0,1)随机数 ；令则第34页，共56页。35第三章 数据添加算法 在Bayes统计或极大似然估计的计算中，经常会遇到这样一类问题：设我们能观测到的数据是Y，关于Y的后验分布p(|Y)很复杂，难以直接进行各种统计计算.假如我们能假定一些没有能观察到的潜在数据Z为已知（譬如，Y为某变量的截尾观测值，Z为该变

15、量的真值），则可能得到一个关于的简单的添加后验分布p(|Y,Z),利用p(|Y,Z)的简单性我们可以进行各种计算，如极大化，抽样等，然后回过头来，又可以对Z的假定做检查或改进。如此进行，我们就将一个复杂的极大化问题转变为一系列简单的极大化或抽样。在统计上，这种处理问题的方法称为“数据添加算法”。 常用的“数据添加算法”有EM算法和Markov Chain Monte Carlo方法。第35页，共56页。363.1 EM算法 先考虑一种简单情形。设某元件的失效时间Y关于变量x有直线回归关系，假设在一次试验中得到一批数据，如图， “”表示该元件失效时间坐标， ”“表示对应元件的截尾时间（小于失效时

16、间）。如果直线斜率和截矩的估计值已知，则我们可以在真实数据不小于截尾数据的前提下将各个被截尾的失效时间估计出来，从而得到所谓的”完全数据“，由此完全数据，重新对直线的斜率及截矩进行估计，再依据新的估计量，得到新的”完全数据“。如此循环往复，则将一个复杂的估计问题替换成一系列简单的估计问题。将之一般化，就给出EM算法。第36页，共56页。37EM算法是一种迭代方法，主要用来求后验分布的众数（即极大似然估计）。它的每一步迭代由两步组成：E步（求期望）和M步（极大化）。一般地，以p(|Y)表示基于Y的的后验密度，称为观测后验分布； p(|Y,Z)表示添加数据Z后得到的的后验密度，称为添加后验分布；

17、p(Z|,Y)表示在给定观测数据Y和参数条件下Z的条件密度。我们的目的是计算p(|Y)的众数。于是EM算法如下进行。记 为第i+1次迭代开始时后验众数的估计值，则第i+1次迭代的两步为E步：将p(|Y,Z)或log p(|Y,Z)关于Z的条件分布求期望，从而把Z积掉，即第37页，共56页。38M步：将 极大化，即找到一个点 ，使将上述E，M步循环进行，直至 充分小为止。例3.1 设总体X的分布律为其中(0,1),现进行了X1234pk197次试验，观察到1,2,3,4的频数为取的先验分布()为U(0,1)分布,则的观察后验分布为第38页，共56页。39现假设X=1可以分解为两部分，其发生概率分

18、别为1/2和/4,令和y1-Z分别表示试验结果中落入这两部分的次数（是不能观测到的潜在数据），则的添加后验分布为(3.1.1)(3.1.2)显然，用(3.1.2)式求极值比(3.1.1)式简单。迭代如下：第39页，共56页。40E步：在给定下，M步：将 关于极大化得可以证明，在关于logp(|Y)的很一般的条件下，由算法得到的估计序列 收敛到的稳定点。（不能保证是极大值点）。较为可行的办法是选几个不同的初值迭代，然后在诸估计值中加以选择，这可减轻初值选取对结果的影响）第40页，共56页。41估计的精度假设EM算法最后的结果是 ，则根据似然估计的渐近正态性，其渐近方差可用 Fisher观测信息的

19、倒数近似。（证明见高等数理统计p126定理2.5.4)第41页，共56页。423.2 Markov Chain Monte Carlo方法对于较简单的后验分布，可直接计算或静态MC等近似计算方法。但在实际中，观测后验分布往往是复杂的，高维的，非标准形式的分布，上述方法都难以实施。对于这类问题，一种简单且行之有效的Bayes计算方法就是MCMC。EM算法得到的是后验分布的众数，有时我们希望得到其它一些后验量如后验均值，方差，后验分布的分位数等。计算这些后验量都可归结为关于后验分布积分的计算。具体地，设 为后验密度，我们要计算的后验量可写成某函数f(x)关于的期望(3.2.1)第42页，共56页。

20、433.2.1 基本思路MCMC方法的基本思想是通过建立一个平稳分布为(x)的Markov链来得到(x)的样本，基于这些样本可以作各种统计推断。比如，若得到了平稳分布为(x)的Markov链的样本轨道 ，则(3.2.1)可估计为(3.2.2)注 由Markov链平稳分布的概念可知，不论Markov链从什么初始状态出发，经过一段时间后，各个时间的边际分布都是平稳分布，因此可将经过某个m时间之后的观察值看作平稳分布(x)的样本。由遍历性定理可知，MCMC的关键是如何构造平稳分布为的Markov链的转移核p(x,y)第43页，共56页。44MCMC方法可概括为如下三步：(1) 在X上选一个“合适”的

21、Markov链，确定其转移核p(x,y),使链的平稳分布为。(2)由X中某一点X(0)出发，用(1)中的Markov链产生序列X1,Xn;(3) 对某个m和大的n，任一函数f(x)的期望估计如下MCMC有许多研究专题，如链的收敛性判断（m大小的确定），链的长度（n的大小）的确定，估计误差等等。以下主要讨论转移核的构造。第44页，共56页。453.2.2 满条件分布MCMC主要用于多变量，非标准形式，且各变量间相互不独立时分布的模拟。令 ，我们总可以写出其中 。如果(3.2.1)式中右端各个因子能够直接模拟，则只需要进行静态模拟（抽样过程中不改变抽样分布）。实际中很难满足上述条件，因此需进行动态

22、模拟（抽样分布随模拟的进行而改变），如MCMC，此时满条件分布扮演了一个重要角色。(3.2.1)在导出满条件分布时，应注意到这样一个事实：记 ，(3.2.2)第45页，共56页。46等价地，若 ，且 ，则(3.2.3)一般地，用y表示观测数据， ，其中 分别表示参数，超参数和缺损数据，则有其中， 表示完全数据的密度函数， 表示先验分布， 表示超参数的分布。有(3.2.2),各变量的满条件分布如下：第46页，共56页。47例3.2.1 设(X1,X2)的联合密度为且 ，则其满条件分布为第47页，共56页。483.2.3 Gibbs 抽样思想：设 的密度为 ，任意固定TN，在给定 条件下，如下定义

23、随机变量 具有密度函数 ，则对任一可测集B，因而X的密度也是 。上述过程定义了一个由X到X的转移核，且其相应的平稳分布是。这样构造的MCMC称为Gibbs抽样。当T只有一个元素时称为单元素Gibbs抽样。第48页，共56页。49单元素Gibbs抽样具体步骤如下：在给定起始点 后，假定第t次迭代开始时的估计值为 ，则第t次迭代分为如下n步：（1）由满条件分布 抽取 (i)由满条件分布 抽取 ; (n)由满条件分布 抽取记 ，则 是平稳分布为的Markov链的实现值，其由x到x的转移概率函数为第49页，共56页。503.2.3 Metroplis-Hastings方法在Gibbs抽样中， 可能很难抽取，这时可采用更一般化的Me