蒙特卡罗方法简介.ppt

《蒙特卡罗方法简介.ppt》由会员分享，可在线阅读，更多相关《蒙特卡罗方法简介.ppt（56页珍藏版）》请在三一文库上搜索。

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

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

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

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

5、x),定义,定理2.1 设随机变量U服从U(0,1)分布，则 的分布函数为F(x). 由定理2.1,要生成分布函数为F(x)的随机数，可先生成U(0,1)随机数U，则可得到随机数X=F-1(U),2.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,1); 2)取 ,确定i，使 3) 由pi(x)抽取x.,2.3 筛选抽样 当

6、p(x)难以直接抽样时，如果可以将p(x) 表示成p(x)=ch(x)g(x),其中h(.)是一密度函数且易于抽样，而0g(y),回到1） 上述方法就是筛选抽样法，它是一种非常重要的抽样方法，可解决许多难以直接抽样的分布的抽样问题。,h(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）如果u p(y)/M(y),回到1）。,筛选抽样的理论依据如下： 定理 设X的密度函数为p(x),且p(x)

7、=ch(x)g(x)，其中0g(x)1,c1 ,h(.)是一密度函数.令U和Y分别服从U(0,1)和h(y)，则在Ug(Y)的条件下，Y的条件密度为,例2.2 设 已知。 注意到,取,则,于是， 的随机数可如下抽取,1）由U(0,1) 抽取u， 2）由h(y)抽取y；（可使用逆变换法） 3）当y(0,1时，如果 ，则x=y， 否则转到1）； 4）当y1时，如果 ，则x=y， 否则转到1）；,2.4 随机向量的抽样法,设X1,Xk的联合概率密度为,定理2.4 设U1,Uk是独立同分布的U(0,1)变量， X1,Xk是方程,的解，其中 是对应于 的分布函数，则X1,Xk的分布为(2.4).,(2.

8、4),(2.5),随机向量的逆变换抽样法： 由U(0,1)分布独立地抽取u1,uk; 用方程（2.5)解x1,xk,例2.3 设X1,X2的联合密度函数为,试生成X1,X2的随机数。,解：,相应的边际分布函数和条件分布函数分别为,方程(2.5)变为,此方程不易解，不妨交换两自变量的次序,相应的边际分布函数和条件分布函数分别为,方程(2.5)变为,对服从特定分布的随机向量有一些特殊的抽样方法。,例2.6 试生成k维正态分布 的随机数。,解：注意到若 ，则存在下三角阵,使,其中C可由迭代实现：,首先，由 ，有,从而,。因,于是,得,依此类推，,一般迭代公式为,至此，我们可以给出k维正态分布的抽样步

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

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

11、f(x)形状接近，应能降低方差。这就是重要抽样法的基本思想。,其方差与g(x)有关。问题变为，如何选择g(.)使估计的方差最小。,例2.5.1 分别用投点法，均值法，重要抽样法，求积分 ,比较各种方法的有效性。,解 i)投点法 1)产生随机数 2) 对每对 ，记 的次数为n0. 则,G,ii)均值法 1)产生随机数 2),iii)重要抽样法 由重要抽样法的思想，需选择一个与 相似的密度函数。由Taylor展开式 取,1)产生随机数 2) 取 则,（数值计算）,模拟结果,二、分层抽样法 另一种利用贡献率大小来降低估计方差的方法是分层抽样法。它首先把样本空间D分成一些不交的小区间 ，然后在各小区间

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

13、易知，在n固定下，当 时,的方差最小，为,如果我们将区间进行10等份，并确定出最优的抽样次数分配： ，则可得到分层抽样法估计的方差为.,一般地，若诸 已知，在n固定下，当 时，估计的方差最小，为,分层抽样法在实施上有两个主要问题，其一是怎样划分区间，简单而常用的方法是将区间等分；另一个问题是在区间划分好后如何确定抽样次数的分配。由于在实际中 总是未知的，因而前面最优分配的结论无法应用。即使如此，分层抽样法还是有其作用的。可以证明，即使取简单的分配 也有,事实上，取 ，代入(2.5.5)得,由Cauchy-Schwarz不等式，有,据此，在(2.5.6)式两端各乘以 并相加得,于是,三、关联抽样

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

15、回过头来，又可以对Z的假定做检查或改进。如此进行，我们就将一个复杂的极大化问题转变为一系列简单的极大化或抽样。在统计上，这种处理问题的方法称为“数据添加算法”。 常用的“数据添加算法”有EM算法和Markov Chain Monte Carlo方法。,3.1 EM算法,先考虑一种简单情形。设某元件的失效时间Y关于变量x有直线回归关系，假设在一次试验中得到一批数据，如图， “”表示该元件失效时间坐标， ”“表示对应元件的截尾时间（小于失效时间）。,如果直线斜率和截矩的估计值已知，则我们可以在真实数据不小于截尾数据的前提下将各个被截尾的失效时间估计出来，从而得到所谓的”完全数据“，由此完全数据，重

16、新对直线的斜率及截矩进行估计，再依据新的估计量，得到新的”完全数据“。如此循环往复，则将一个复杂的估计问题替换成一系列简单的估计问题。将之一般化，就给出EM算法。,EM算法是一种迭代方法，主要用来求后验分布的众数（即极大似然估计）。它的每一步迭代由两步组成：E步（求期望）和M步（极大化）。 一般地，以p(|Y)表示基于Y的的后验密度，称为观测后验分布； p(|Y,Z)表示添加数据Z后得到的的后验密度，称为添加后验分布； p(Z|,Y)表示在给定观测数据Y和参数条件下Z的条件密度。 我们的目的是计算p(|Y)的众数。 于是EM算法如下进行。记 为第i+1次迭代开始时后验众数的估计值，则第i+1次

17、迭代的两步为 E步：将p(|Y,Z)或log p(|Y,Z)关于Z的条件分布求期望，从而把Z积掉，即,M步：将 极大化，即找到一个点 ，使,将上述E，M步循环进行，直至 充分小为止。,例3.1 设总体X的分布律为 其中(0,1),现进行了,197次试验，观察到1,2,3,4的频数为 取的先验分布()为U(0,1)分布,则的观察后验分布为,现假设X=1可以分解为两部分，其发生概率分别为1/2和/4,令和y1-Z分别表示试验结果中落入这两部分的次数（是不能观测到的潜在数据），则的添加后验分布为,(3.1.1),(3.1.2),显然，用(3.1.2)式求极值比(3.1.1)式简单。迭代如下：,E步：

18、在给定下，,M步：将 关于极大化得,可以证明，在关于logp(|Y)的很一般的条件下，由算法得到的估计序列 收敛到的稳定点。（不能保证是极大值点）。较为可行的办法是选几个不同的初值迭代，然后在诸估计值中加以选择，这可减轻初值选取对结果的影响）,估计的精度,假设EM算法最后的结果是 ，则根据似然估计的渐近正态性，其渐近方差可用 Fisher观测信息的倒数,近似。（证明见高等数理统计p126定理2.5.4),3.2 Markov Chain Monte Carlo方法,对于较简单的后验分布，可直接计算或静态MC等近似计算方法。但在实际中，观测后验分布往往是复杂的，高维的，非标准形式的分布，上述方法

19、都难以实施。对于这类问题，一种简单且行之有效的Bayes计算方法就是MCMC。,EM算法得到的是后验分布的众数，有时我们希望得到其它一些后验量如后验均值，方差，后验分布的分位数等。计算这些后验量都可归结为关于后验分布积分的计算。具体地，设 为后验密度，我们要计算的后验量可写成某函数f(x)关于的期望,(3.2.1),3.2.1 基本思路,MCMC方法的基本思想是通过建立一个平稳分布为(x)的Markov链来得到(x)的样本，基于这些样本可以作各种统计推断。比如，若得到了平稳分布为(x)的Markov链的样本轨道 ，则(3.2.1)可估计为,(3.2.2),注 由Markov链平稳分布的概念可知

20、，不论Markov链从什么初始状态出发，经过一段时间后，各个时间的边际分布都是平稳分布，因此可将经过某个m时间之后的观察值看作平稳分布(x)的样本。由遍历性定理可知，,MCMC的关键是如何构造平稳分布为的Markov链的转移核p(x,y),MCMC方法可概括为如下三步： (1) 在X上选一个“合适”的Markov链，确定其转移核p(x,y),使链的平稳分布为。 (2)由X中某一点X(0)出发，用(1)中的Markov链产生序列X1,Xn; (3) 对某个m和大的n，任一函数f(x)的期望估计如下,MCMC有许多研究专题，如链的收敛性判断（m大小的确定），链的长度（n的大小）的确定，估计误差等等

21、。以下主要讨论转移核的构造。,3.2.2 满条件分布,MCMC主要用于多变量，非标准形式，且各变量间相互不独立时分布的模拟。,令 ，我们总可以写出,其中 。如果(3.2.1)式中右端各个因子能够直接模拟，则只需要进行静态模拟（抽样过程中不改变抽样分布）。实际中很难满足上述条件，因此需进行动态模拟（抽样分布随模拟的进行而改变），如MCMC，此时满条件分布扮演了一个重要角色。,(3.2.1),在导出满条件分布时，应注意到这样一个事实： 记 ，,(3.2.2),等价地，若 ，且 ，则,(3.2.3),一般地，用y表示观测数据， ，其中 分别表示参数，超参数和缺损数据，则有,其中， 表示完全数据的密度

22、函数， 表示先验分布， 表示超参数的分布。有(3.2.2),各变量的满条件分布如下：,例3.2.1 设(X1,X2)的联合密度为,且 ，则其满条件分布为,3.2.3 Gibbs 抽样,思想：设 的密度为 ，任意固定TN，在给定 条件下，如下定义随机变量 具有密度函数 ，则对任一可测集B，,因而X的密度也是 。,上述过程定义了一个由X到X的转移核，且其相应的平稳分布是。这样构造的MCMC称为Gibbs抽样。当T只有一个元素时称为单元素Gibbs抽样。,单元素Gibbs抽样具体步骤如下： 在给定起始点 后，假定第t次迭代开始时的估计值为 ，则第t次迭代分为如下n步： （1）由满条件分布 抽取 (i

23、)由满条件分布 抽取 ; (n)由满条件分布 抽取,记 ，则 是平稳分布为的Markov链的实现值，其由x到x的转移概率函数为,3.2.3 Metroplis-Hastings方法,在Gibbs抽样中， 可能很难抽取，这时可采用更一般化的Metroplis-Hastings方法。其思路如下： 任意选择一个不可约的转移概率q(.,.)以及一个函数(.,.)，01,对任一组合(x,x)(xx),定义,则p(x,x)形成一个转移核。,此方法的实施比较直观，首先由q(x,x)产生一个潜在的转移x x，然后根据概率(x,x)决定是否转移。于是，在有了x之后，可产生一U(0,1)随机数u,如果u(x,x)

24、,则X(i+1)=x,否则， X(i+1)=x。,注意：在有了q(.,.)后，应选择(.,.)使相应的p(x,x)以(x)为平稳分布。 定理3.2.1 给定q(x,x),取,则由以下p(x,x)产生的Markov链是可逆的，且以(x)为平稳分布,由定理3.2.1可以看出，建议分布q(x,x)可取各种形式。 常用的建议分布选择方法： (1)Metropolis选择：对称的建议分布，即,此时,例如,（2）独立抽样,取,则,一般，要使独立抽样有好的效果，q(x)应接近(x),比较安全的办法是使q(x)的尾比(x)重。,(3)单元素Metropolis-Hastings算法,产生向量X的样本有时是困难的，通常是利用满条件分布对X的分量进行逐个抽样。,考虑 的条件分布，选择一个转移核 ，固定 不变，由q 产生一个可能的 ，然后以概率决定是否接受其作为下一状态。,（续）例3.2.1 设(X1,X2)的联合密度为,试产生 的后验分布样本。,解：由例3.2.1, 各分量的满条件分布为,取,则抽样步骤为： 选初值 按 产生随机数 按 U(0,1)产生随机数 产生一U(0,1)随机数u,如果 则否则 ，否则,