蒙特卡罗方法课件

《蒙特卡罗方法课件》由会员分享，可在线阅读，更多相关《蒙特卡罗方法课件（56页珍藏版）》请在装配图网上搜索。

1、目目 录录 第一章蒙特卡罗方法概述第二章随机数的产生第三章EM算法和MCMC方法参考书参考书:1.1.茆诗松等茆诗松等,高等数理统计高等数理统计(第第6 6章章),),高等教育出版社高等教育出版社,1998;,1998;2.2.徐钟济，蒙特卡罗方法，上海科学技术出版社徐钟济，蒙特卡罗方法，上海科学技术出版社1PPT学习交流第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述 蒙特卡罗方法又称随机抽样技巧或统计试蒙特卡罗方法又称随机抽样技巧或统计试验方法。验方法。蒙特卡罗方法是一种计算方法，但与一般蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它以概率统计理论数值计算方法有很大区别。它以

2、概率统计理论为基础。由于蒙特卡罗方法能够比较逼真地描为基础。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域方法难以解决的问题，因而该方法的应用领域日趋广泛。日趋广泛。2PPT学习交流1.蒙特卡罗方法的基本思想蒙特卡罗方法的基本思想 理论基础：大数定律；中心极限定理；F(X)U(0,1)。基本思想：1.当所求问题的解是某个事件的概率，或者是某个随机变量的期望，或与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或该随机变量若干个观察值的算术平均值，根据大数定律得到问题的解

3、；2.要生成分布函数为F(x)的随机数，可先生成U(0,1)随机数F，则可得到随机数X=F-1(F)。3PPT学习交流例（利用例（利用MC进行欧式期权定价）设股票价格进行欧式期权定价）设股票价格St服服从风险中性测度下的几何从风险中性测度下的几何Brown运动：运动：其离散化形式为其离散化形式为根据金融工程理论，设现在股票价格为根据金融工程理论，设现在股票价格为S0，T时时刻到期（单位天），敲定价为刻到期（单位天），敲定价为K的欧式看涨期权的欧式看涨期权的价格为的价格为MC方案：按照（方案：按照（1）递推产生）递推产生n条风险中性测度下的条风险中性测度下的轨道，提取出轨道，提取出ST(n)；（

4、；（2）4PPT学习交流2.蒙特卡罗方法的蒙特卡罗方法的误差 根据中心极限定理如果随机变量序列根据中心极限定理如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差独立同分布，且具有有限非零的方差2，即，即则则当N充分大时，有如下的近似式它表明，误差收敛速度的阶为它表明，误差收敛速度的阶为 以概率以概率1-成立。成立。5PPT学习交流通常，蒙特卡罗方法的误差通常，蒙特卡罗方法的误差定义为定义为关于蒙特卡罗方法的误差需说明两点：关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。计算方法是

5、有区别的。第二，误差中的均方差第二，误差中的均方差是未知的，必须使用其估计值是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出来代替，在计算所求量的同时，可计算出 。6PPT学习交流减小方差的各种技巧 显然，当给定置信度显然，当给定置信度后，误差后，误差由由和和N N决定。决定。要减小要减小，或者是增大，或者是增大N N，或者是减小方差，或者是减小方差2 2。在。在固定的情况下，要把精度提高一个数量级，试验次数固定的情况下，要把精度提高一个数量级，试验次数N N需增加两个数量级。因此，单纯增大需增加两个数量级。因此，单纯增大N N不是一个有效不是一个有效的办法。降低方差的各种技巧，

6、引起了人们的普遍注的办法。降低方差的各种技巧，引起了人们的普遍注意。意。一般来说，降低方差的技巧，往往会使观察一个子样一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中特卡罗方法中效率效率的概念。的概念。它定义为 其中其中c c是观察一个子样的平均费用。是观察一个子样的平均费用。7PPT学习交流3.蒙特卡罗方法的特

7、点蒙特卡罗方法的特点优点优点1)1)能够比较逼真地描述具有能够比较逼真地描述具有随机性质的事物的特点及随机性质的事物的特点及物理实验过程。物理实验过程。2)2)受几何条件限制小。受几何条件限制小。3)3)收敛速度与问题的维数无收敛速度与问题的维数无关。关。4)4)误差容易确定。误差容易确定。5)5)程序结构简单，易于实现。程序结构简单，易于实现。缺点缺点1)收敛速度慢。收敛速度慢。2)误差具有概率性。误差具有概率性。8PPT学习交流第二章随机数的产生2.1 2.1 逆变换法逆变换法设随机变量设随机变量X X的分布函数为的分布函数为F(x),F(x),定义定义定理定理2.1 设随机变量设随机变量

8、U服从服从U(0,1)分布，则分布，则的分布函数为的分布函数为F(x).由定理由定理2.1,要生成分布函数为要生成分布函数为F(x)的随机数，可先生的随机数，可先生成成U(0,1)随机数随机数U，则可得到随机数，则可得到随机数X=F-1(U)9PPT学习交流2.2 合成法合成法如果如果X的密度函数的密度函数p(x)难于抽样难于抽样,而而X关于关于Y的条件的条件密度函数密度函数p(x|y)以及以及Y的密度函数的密度函数g(y)均易于抽样，均易于抽样，则则X 的随机数可如下产生：的随机数可如下产生：Step1 由由Y的分布的分布g(y)抽取抽取y;Step2 由由X关于关于Y的条件密度函数的条件密

9、度函数p(x|y)抽取抽取x.例例2.1 设设X的密度函数为的密度函数为由合成法，由合成法，X的随机数可如下抽取：的随机数可如下抽取：1）取）取uU(0,1);2)取取 ,确定确定i，使，使3)由由pi(x)抽取抽取x.10PPT学习交流2.3 筛选抽样筛选抽样 当当p(x)难以直接抽样时，如果可以将难以直接抽样时，如果可以将p(x)表示成表示成p(x)=c h(x)g(x),其中其中h(.)是一密度函数且易于抽样，是一密度函数且易于抽样，而而0g(y),回到回到1）上述方法就是筛选抽样法，它是一种非常重要的抽样上述方法就是筛选抽样法，它是一种非常重要的抽样方法，可解决许多难以直接抽样的分布的

10、抽样问题。方法，可解决许多难以直接抽样的分布的抽样问题。11PPT学习交流h(x)的的选取有多种方法。一种直观的方法是：如果的的选取有多种方法。一种直观的方法是：如果存在一个函数存在一个函数M(x)，满足，满足p(x)M(x),且且令令h(x)=M(x)/c,若若h(x)易于抽样，则易于抽样，则筛选抽样变为筛选抽样变为1）由）由U(0,1)抽取抽取u，由，由h(y)抽取抽取y；2）如果）如果u p(y)/M(y),则则x=y停止；停止；3）如果）如果u p(y)/M(y),回到回到1）。）。筛选抽样的理论依据如下：筛选抽样的理论依据如下：定理定理 设设X的密度函数为的密度函数为p(x),且且p

11、(x)=c h(x)g(x)，其，其中中01时，如果时，如果 ，则，则x=y，否则转到否则转到1）；）；14PPT学习交流2.4 随机向量的抽样法随机向量的抽样法设设X1,Xk的联合概率密度为的联合概率密度为定理定理2.4 设设U1,Uk是独立同分布的是独立同分布的U(0,1)变量，变量，X1,Xk是方程是方程的解，其中的解，其中 是对应于是对应于 的分布函数，则的分布函数，则X1,Xk的分布为的分布为(2.4).(2.4)(2.5)15PPT学习交流随机向量的逆变换抽样法：随机向量的逆变换抽样法：1)由由U(0,1)分布独立地抽取分布独立地抽取u1,uk;2)用方程（用方程（2.5)解解x1

12、,xk例例2.3 设设X1,X2的联合密度函数为的联合密度函数为试生成试生成X1,X2的随机数。的随机数。解：解：16PPT学习交流相应的边际分布函数和条件分布函数分别为相应的边际分布函数和条件分布函数分别为方程方程(2.5)变为变为此方程不易解，不妨交换两此方程不易解，不妨交换两自变量的次序自变量的次序17PPT学习交流相应的边际分布函数和条件分布函数分别为相应的边际分布函数和条件分布函数分别为方程方程(2.5)变为变为对服从特定分布的随机向量有一些特殊的抽样方法。对服从特定分布的随机向量有一些特殊的抽样方法。18PPT学习交流例例2.6 试生成试生成k维正态分布维正态分布 的随机数。的随机

13、数。解：注意到若解：注意到若 ，则存在，则存在下三角阵下三角阵使使其中其中C可由迭代实现：可由迭代实现：首先，由首先，由 ，有，有从而从而。因。因于是于是得得依此类推，依此类推，19PPT学习交流一般迭代公式为一般迭代公式为至此，我们可以给出至此，我们可以给出k维正态分布的抽样步骤：维正态分布的抽样步骤：1)迭代计算迭代计算 ；2)由由N(0,1)分布独立抽取分布独立抽取k个随机数个随机数 ；3)计算计算20PPT学习交流2.5 随机模拟计算随机模拟计算2.5.1 随机投点法随机投点法考虑积分考虑积分 ，设，设a,b有限，有限，0 f(x)M,令令=(x,y):a=(x,y):a x x b,

14、0b,0 y y M,M,并设并设(X,Y)(X,Y)是在是在上均上均匀分布的二维随机向量，其联合密度函数为匀分布的二维随机向量，其联合密度函数为则易见，则易见，是是中曲线中曲线f(x)f(x)下方面积下方面积。假设我们向假设我们向中投点，若点落在中投点，若点落在y=f(x)下方称为中的，下方称为中的，则点中的概率为则点中的概率为21PPT学习交流若我们进行了若我们进行了n次投点，其中次投点，其中n0次中的，则可以得次中的，则可以得到到一个估计一个估计不难看出，不难看出，是是的无偏估计，且其方差为的无偏估计，且其方差为(2.5.1)22PPT学习交流2.5.2 样本均值法样本均值法于是，积分于

15、是，积分注意到，若注意到，若XU(a,b),则则由大数定律，若由大数定律，若 ，则，则MC方法为：方法为：1)独立产生独立产生n个个U(a,b)随机数随机数2)按按(2.5.2)估计估计。(2.5.2)23PPT学习交流可证，在可证，在0 f(x)M条件下，条件下，2.5.3 降低方差的技术降低方差的技术Monte Carlo 方法中一类重要的研究课题是考虑一方法中一类重要的研究课题是考虑一些降低估计方差的技术。常用的方法有：重要抽样些降低估计方差的技术。常用的方法有：重要抽样法，分层抽样法，关联抽样法等。法，分层抽样法，关联抽样法等。一一 重要抽样法重要抽样法由上节，样本平均法比投点法有效，

16、将样本平均法做由上节，样本平均法比投点法有效，将样本平均法做更一般的推广，设更一般的推广，设g(x)是是(a,b)上的密度函数，改写上的密度函数，改写24PPT学习交流由大数定律，若由大数定律，若 ，则，则MC方法为：方法为：1)选择适当的选择适当的g(x),独立产生独立产生n个个g(x)随机数随机数2)由由(2.5.3)估计估计。显然显然(2.5.3)25PPT学习交流从理论上看，因，若f(x)0,取则有因为未知，这是作不到的，但它提示我们取g(x)与f(x)形状接近，应能降低方差。这就是重要抽样法的基本思想。其方差与其方差与g(x)有关。问题变为，如何选择有关。问题变为，如何选择g(.)使

17、估使估计的方差最小。计的方差最小。26PPT学习交流例例2.5.1 分别用投点法，均值法，重要抽样法，求积分别用投点法，均值法，重要抽样法，求积分分 ,比较各种方法的有效性。比较各种方法的有效性。解解 i)投点法投点法1)产生随机数 2)对每对 ，记 的次数为n0.则Gii)均值法均值法1)产生随机数 2)27PPT学习交流iii)重要抽样法重要抽样法由重要抽样法的思想，需选择一个与 相似的密度函数。由Taylor展开式 取1)产生随机数 2)取则（数值计算）真值真值投点法投点法均值法均值法重要抽样法重要抽样法1.718281.87561.81451.7219模拟结果28PPT学习交流二、分层

18、抽样法二、分层抽样法另一种利用贡献率大小来降低估计方差的方法是分层另一种利用贡献率大小来降低估计方差的方法是分层抽样法。它首先把样本空间抽样法。它首先把样本空间D分成一些不交的小区间分成一些不交的小区间 ，然后在各小区间内的抽样数由其贡献大小，然后在各小区间内的抽样数由其贡献大小决定。即，定义决定。即，定义 ，则，则Di内的抽样内的抽样数数ni应与应与pi成正比。成正比。考虑积分考虑积分将将0,1分成m个小区间：则则记记 为第为第i个小区间的长度，个小区间的长度，i=1,m.在每个小区间上的积分值可用均值法估计出来，然在每个小区间上的积分值可用均值法估计出来，然后将其相加即可给出后将其相加即可

19、给出的一个估计。具体步骤为：的一个估计。具体步骤为：29PPT学习交流1)独立产生独立产生U(0,1)随机数随机数2)计算计算3)计算计算于是可得于是可得的估计为的估计为(2.5.4)易见，易见，是是的无偏估计，其方差为的无偏估计，其方差为(2.5.5)(2.5.6)30PPT学习交流续例续例2.5.1 考察分层抽样法求积分考察分层抽样法求积分 的方的方差。差。解：先将区间解：先将区间0,1划分成两个小区间划分成两个小区间0,0.5,0.5,1,则则设一共抽设一共抽n个随机数，其中在个随机数，其中在0,0.5)上抽上抽n1个，则个，则使用分层抽样法求得使用分层抽样法求得 的方差为的方差为31P

20、PT学习交流对对n1求导易知，在求导易知，在n固定下，当固定下，当 时时的方差最小，为的方差最小，为如果我们将区间进行如果我们将区间进行10等份，并确定出最优的抽样等份，并确定出最优的抽样次数分配：次数分配：，则可得到，则可得到分层抽样法估计的方差为分层抽样法估计的方差为.一般地，若诸一般地，若诸 已知，在已知，在n固定下，当固定下，当时，估计的方差最小，为时，估计的方差最小，为32PPT学习交流分层抽样法在实施上有两个主要问题，其一是怎样划分区间，简单而常用的方法是将区间等分；另一个问题是在区间划分好后如何确定抽样次数的分配。由于在实际中 总是未知的，因而前面最优分配的结论无法应用。即使如此

21、，分层抽样法还是有其作用的。可以证明，即使取简单的分配也有事实上，取 ，代入(2.5.5)得由Cauchy-Schwarz不等式，有据此，在(2.5.6)式两端各乘以 并相加得 于是于是33PPT学习交流三、关联抽样法三、关联抽样法考虑积分差考虑积分差若用若用 估计估计，则其方差为，则其方差为显然，在显然，在 确定后，确定后，正相关度越高，则正相关度越高，则 的方差越小。这便是关联抽样的方差越小。这便是关联抽样法的基本出发点。法的基本出发点。考虑用重要抽样法来估计考虑用重要抽样法来估计I1,I2,即改写即改写为为1)产生产生n个个U(0,1)随机数随机数 ；令；令2)则则34PPT学习交流第三

22、章数据添加算法 在Bayes统计或极大似然估计的计算中，经常会遇到这样一类问题：设我们能观测到的数据是Y，关于Y的后验分布p(|Y)很复杂，难以直接进行各种统计计算.假如我们能假定一些没有能观察到的潜在数据Z为已知（譬如，Y为某变量的截尾观测值，Z为该变量的真值），则可能得到一个关于的简单的添加后验分布p(|Y,Z),利用p(|Y,Z)的简单性我们可以进行各种计算，如极大化，抽样等，然后回过头来，又可以对Z的假定做检查或改进。如此进行，我们就将一个复杂的极大化问题转变为一系列简单的极大化或抽样。在统计上，这种处理问题的方法称为“数据添加算法”。常用的“数据添加算法”有EM算法和Markov C

23、hain Monte Carlo方法。35PPT学习交流3.1EM算法 先考虑一种简单情形。设某元件的失效时间Y关于变量x有直线回归关系，假设在一次试验中得到一批数据，如图，“”表示该元件失效时间坐标，”“表示对应元件的截尾时间（小于失效时间）。如果直线斜率和截矩的估计值已知，则我们可以在真实数据不小于截尾数据的前提下将各个被截尾的失效时间估计出来，从而得到所谓的”完全数据“，由此完全数据，重新对直线的斜率及截矩进行估计，再依据新的估计量，得到新的”完全数据“。如此循环往复，则将一个复杂的估计问题替换成一系列简单的估计问题。将之一般化，就给出EM算法。36PPT学习交流EM算法是一种迭代方法，

24、主要用来求后验分布的众算法是一种迭代方法，主要用来求后验分布的众数（即极大似然估计）。它的每一步迭代由两步组成：数（即极大似然估计）。它的每一步迭代由两步组成：E步（求期望）和步（求期望）和M步（极大化）。步（极大化）。一般地，以一般地，以p(|Y)|Y)表示表示基于基于Y Y的的后验密度，称为的的后验密度，称为观测后验分布；观测后验分布；p(|Y,Z)表示添加数据表示添加数据Z后得到的后得到的的后验密度，称为添加后验分布；的后验密度，称为添加后验分布；p(Z|,Y)表示在给表示在给定观测数据定观测数据Y和参数和参数条件下条件下Z Z的条件密度。的条件密度。我们的目的是计算我们的目的是计算p(

25、|Y)的众数。的众数。于是于是EM算法如下进行。记算法如下进行。记 为第为第i+1次迭代开始时次迭代开始时后验众数的估计值，则第后验众数的估计值，则第i+1次迭代的两步为次迭代的两步为E步：将步：将p(|Y,Z)或或log p(|Y,Z)关于关于Z的条件分布求期的条件分布求期望，从而把望，从而把Z积掉，即积掉，即37PPT学习交流M步：将步：将 极大化，即找到一个点极大化，即找到一个点 ，使，使将上述将上述E，M步循环进行，直至步循环进行，直至 充分小为充分小为止。止。例例3.1 设总体设总体X的分布律为的分布律为其中其中(0,1),(0,1),现进行了现进行了X1234pk197次试验，观察

26、到次试验，观察到1,2,3,4的频数为的频数为取取的先验分布的先验分布()为为U(0,1)U(0,1)分布分布,则则的观察后的观察后验分布为验分布为38PPT学习交流现假设现假设X=1可以分解为两部分，其发生概率分别为可以分解为两部分，其发生概率分别为1/2和和/4,/4,令和令和y y1 1-Z-Z分别表示试验结果中落入这分别表示试验结果中落入这两部分的次数（是不能观测到的潜在数据），则两部分的次数（是不能观测到的潜在数据），则的添加后验分布为的添加后验分布为(3.1.1)(3.1.2)显然，用显然，用(3.1.2)式求极值比式求极值比(3.1.1)式简单。迭代如下：式简单。迭代如下：39P

27、PT学习交流E步：在给定下，步：在给定下，M步：将步：将 关于关于极大化得极大化得可以证明，在关于logp(|Y)的很一般的条件下，由算法得到的估计序列 收敛到的稳定点。（不能保证是极大值点）。较为可行的办法是选几个不同的初值迭代，然后在诸估计值中加以选择，这可减轻初值选取对结果的影响）40PPT学习交流估计的精度假设假设EM算法最后的结果是算法最后的结果是 ，则根据似然估计的渐，则根据似然估计的渐近正态性，其渐近方差可用近正态性，其渐近方差可用 Fisher观测信息的倒数观测信息的倒数近似。近似。（证明见高等数理统计p126定理2.5.4)41PPT学习交流3.2MarkovChainMon

28、teCarlo方法对于较简单的后验分布，可直接计算或静态对于较简单的后验分布，可直接计算或静态MC等近等近似计算方法。但在实际中，观测后验分布往往是复杂似计算方法。但在实际中，观测后验分布往往是复杂的，高维的，非标准形式的分布，上述方法都难以实的，高维的，非标准形式的分布，上述方法都难以实施。对于这类问题，一种简单且行之有效的施。对于这类问题，一种简单且行之有效的Bayes计计算方法就是算方法就是MCMC。EM算法得到的是后验分布的众数，有时我们希望得算法得到的是后验分布的众数，有时我们希望得到其它一些后验量如后验均值，方差，后验分布的到其它一些后验量如后验均值，方差，后验分布的分位数等。计算

29、这些后验量都可归结为关于后验分分位数等。计算这些后验量都可归结为关于后验分布积分的计算。具体地，设布积分的计算。具体地，设 为后验密度，为后验密度，我们要计算的后验量可写成某函数我们要计算的后验量可写成某函数f(x)关于关于的期望的期望(3.2.1)42PPT学习交流3.2.1基本思路MCMC方法的基本思想是通过建立一个平稳分布为方法的基本思想是通过建立一个平稳分布为(x)(x)的的MarkovMarkov链来得到链来得到(x)(x)的样本，基于这些样本的样本，基于这些样本可以作各种统计推断。比如，若得到了平稳分布为可以作各种统计推断。比如，若得到了平稳分布为(x)(x)的的Markov链的链

30、的样本轨道样本轨道 ，则，则(3.2.1)(3.2.1)可可估计为估计为(3.2.2)注注 由Markov链平稳分布的概念可知，不论Markov链从什么初始状态出发，经过一段时间后，各个时间的边际分布都是平稳分布，因此可将经过某个m时间之后的观察值看作平稳分布(x)的样本。由遍历性定理可知，MCMC的关键是如何构造平稳分布为的关键是如何构造平稳分布为的的Markov链的转移核p(x,y)43PPT学习交流MCMC方法可概括为如下三步：方法可概括为如下三步：(1)在在X X上选一个上选一个“合适合适”的的Markov链，确定其转移链，确定其转移核核p(x,y),使链的平稳分布为使链的平稳分布为。

31、(2)(2)由由X X中某一点中某一点X X(0)(0)出发，用出发，用(1)(1)中的中的Markov链产链产生序列生序列X1,Xn;(3)对某个对某个m和大的和大的n，任一函数，任一函数f(x)的期望估计如下的期望估计如下MCMC有许多研究专题，如链的收敛性判断（m大小的确定），链的长度（n的大小）的确定，估计误差等等。以下主要讨论转移核的构造。44PPT学习交流3.2.2满条件分布MCMC主要用于多变量，非标准形式，且各变量间主要用于多变量，非标准形式，且各变量间相互不独立时分布的模拟。相互不独立时分布的模拟。令令 ，我们总可以写出，我们总可以写出其中 。如果(3.2.1)式中右端各个因

32、子能够直接模拟，则只需要进行静态模拟（抽样过程中不改变抽样分布）。实际中很难满足上述条件，因此需进行动态模拟（抽样分布随模拟的进行而改变），如MCMC，此时满条件分布扮演了一个重要角色。(3.2.1)在导出满条件分布时，应注意到这样一个事实：记 ，(3.2.2)45PPT学习交流等价地，若等价地，若 ，且，且 ，则，则(3.2.3)一般地，用一般地，用y表示观测数据，表示观测数据，其中，其中 分分别表示参数，超参数和缺损数据，则有别表示参数，超参数和缺损数据，则有其中，其中，表示完全数据的密度函数，表示完全数据的密度函数，表表示先验分布，示先验分布，表示超参数的分布。有表示超参数的分布。有(3

33、.2.2),各变各变量的满条件分布如下：量的满条件分布如下：46PPT学习交流例例3.2.1 设设(X1,X2)的联合密度为的联合密度为且且 ，则其满条件分布为，则其满条件分布为47PPT学习交流3.2.3Gibbs抽样思想：设思想：设 的密度为的密度为 ，任意，任意固定固定T N，在给定，在给定 条件下，如下定义随条件下，如下定义随机变量机变量 具有具有密度函数密度函数 ，则对任一可测集，则对任一可测集B，因而因而X的密度也是的密度也是 。上述过程定义了一个由上述过程定义了一个由X到到X的转移核，且其相应的转移核，且其相应的平稳分布是的平稳分布是。这样构造的。这样构造的MCMCMCMC称为称

34、为GibbsGibbs抽样。抽样。当当T T只有一个元素时称为单元素只有一个元素时称为单元素Gibbs抽样。抽样。48PPT学习交流单元素单元素Gibbs抽样具体步骤如下：抽样具体步骤如下：在给定起始点在给定起始点 后，假定第后，假定第t次迭代次迭代开始时的估计值为开始时的估计值为 ，则第，则第t次迭代分为如下次迭代分为如下n步：步：（1）由满条件分布）由满条件分布 抽取抽取(i)由满条件分布由满条件分布 抽取抽取 ;(n)由满条件分布由满条件分布 抽取抽取记记 ，则，则 是平稳分布为是平稳分布为的的Markov链的实现值，其由链的实现值，其由x到到x的转移概率函数为的转移概率函数为49PPT

35、学习交流3.2.3Metroplis-Hastings方法在在Gibbs抽样中，抽样中，可能很难抽取，这时可采用可能很难抽取，这时可采用更一般化的更一般化的Metroplis-Hastings方法。其思路如下：方法。其思路如下：任意选择一个不可约的转移概率任意选择一个不可约的转移概率q(.,.)以及一个以及一个函数函数(.,.)(.,.)，00 1,1,对任一组合对任一组合(x,x(x,x)(x)(x x x),),定义定义则则p(x,x)形成一个转移核。形成一个转移核。此方法的实施比较直观，首先由此方法的实施比较直观，首先由q(x,x)产生一个潜在产生一个潜在的转移的转移x x，然后根据概率

36、，然后根据概率(x,x)决定是否转移。于决定是否转移。于是，在有了是，在有了x x之后，可产生一之后，可产生一U(0,1)U(0,1)随机数随机数u,u,如果如果u u(x,x(x,x),),则则X X(i+1)(i+1)=x=x,否则，否则，X(i+1)=x。50PPT学习交流注意：在有了注意：在有了q(.,.)q(.,.)后，应选择后，应选择(.,.)(.,.)使相应的使相应的p(x,xp(x,x)以以(x)为平稳分布。为平稳分布。定理定理3.2.1 3.2.1 给定给定q(x,xq(x,x),),取取则由以下则由以下p(x,x)产生的产生的Markov链是可逆的，且以链是可逆的，且以(x

37、)为平稳分布为平稳分布51PPT学习交流由定理由定理3.2.1可以看出，建议分布可以看出，建议分布q(x,x)可取各种形式。可取各种形式。常用的建议分布选择方法：常用的建议分布选择方法：(1)Metropolis选择：对称的建议分布，即选择：对称的建议分布，即此时此时例如例如52PPT学习交流（2）独立抽样）独立抽样取取则则一般，要使独立抽样有好的效果，一般，要使独立抽样有好的效果，q(x)应接近应接近(x),(x),比较安全的办法是使比较安全的办法是使q(x)q(x)的尾比的尾比(x)重。重。(3)单元素单元素Metropolis-Hastings算法算法产生向量产生向量X的样本有时是困难的

38、，通常是利用满条的样本有时是困难的，通常是利用满条件分布对件分布对X的分量进行逐个抽样。的分量进行逐个抽样。考虑考虑 的条件分布，选择一个转移核的条件分布，选择一个转移核 ，固定，固定 不变，由不变，由q 产生产生一个可能的一个可能的 ，然后以概率，然后以概率决定是否接受其作决定是否接受其作为下一状态。为下一状态。53PPT学习交流（续）例（续）例3.2.1 设设(X1,X2)的联合密度为的联合密度为试产生试产生 的后验分布样本。的后验分布样本。解：由例解：由例3.2.1,各分量的满条件分布为各分量的满条件分布为54PPT学习交流取取则抽样步骤为：则抽样步骤为：(1)选初值选初值(2)按按 产生随机数产生随机数(3)按按 U(0,1)产生随机数产生随机数(4)(4)产生一产生一U(0,1)U(0,1)随机数随机数u,u,如果如果 则否则则否则 ，否则，否则55PPT学习交流此课件下载可自行编辑修改，供参考！此课件下载可自行编辑修改，供参考！感谢您的支持，我们努力做得更好！感谢您的支持，我们努力做得更好！56PPT学习交流