1、1目目 录录 第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述 第二章第二章 随机数的产生随机数的产生 第三章第三章 EM算法和算法和MCMC方法方法参考书参考书 : 1.1. 茆诗松等茆诗松等, , 高等数理统计高等数理统计( (第第6 6章章), ), 高等教育出版社高等教育出版社,1998;,1998;2.2.徐钟济，蒙特卡罗方法，上海科学技术出版社徐钟济，蒙特卡罗方法，上海科学技术出版社2第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述 蒙特卡罗方法又称随机抽样技巧或统计试蒙特卡罗方法又称随机抽样技巧或统计试验方法。验方法。 蒙特卡罗方法是一种计算方法，但与一般蒙特卡罗方法是一种计算方法，
2、但与一般数值计算方法有很大区别。它以概率统计理论数值计算方法有很大区别。它以概率统计理论为基础。由于蒙特卡罗方法能够比较逼真地描为基础。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域方法难以解决的问题，因而该方法的应用领域日趋广泛。日趋广泛。31.蒙特卡罗方法的基本思想蒙特卡罗方法的基本思想 理论基础：大数定律；中心极限定理； F(X)U(0,1)。基本思想：1.当所求问题的解是某个事件的概率，或者是某个随机变量的期望，或与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率
3、，或该随机变量若干个观察值的算术平均值，根据大数定律得到问题的解；2. 要生成分布函数为F(x)的随机数，可先生成U(0,1)随机数F，则可得到随机数X=F-1(F) 。4例（利用例（利用MC进行欧式期权定价）设股票价格进行欧式期权定价）设股票价格St服服从风险中性测度下的几何从风险中性测度下的几何Brown运动：运动：ttttdSrS dtS dB其离散化形式为其离散化形式为1(1)(0,1)(1)iiiiiSr SSBBN根据金融工程理论，设现在股票价格为根据金融工程理论，设现在股票价格为S0，T时时刻到期（单位天），敲定价为刻到期（单位天），敲定价为K的欧式看涨期权的欧式看涨期权的价格为
4、的价格为rTTCeESKMC方案：按照（方案：按照（1）递推产生）递推产生n条风险中性测度下的条风险中性测度下的轨道，提取出轨道，提取出ST (n)；（；（2） 11nirTTiCeSKn52. 蒙特卡罗方法的蒙特卡罗方法的误差 2/202( )12ztNzP XE XedtN dtexXEXNPxxtNN2/221)(limdxxfXEx)()(022根据中心极限定理如果随机变量序列根据中心极限定理如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差独立同分布，且具有有限非零的方差2 ，即，即则则当N充分大时，有如下的近似式它表明，误差收敛速度的阶为它表明，误差收敛速度的阶为 以
5、概率以概率1-成立成立。6zN通常，蒙特卡罗方法的误差通常，蒙特卡罗方法的误差定义为定义为关于蒙特卡罗方法的误差需说明两点：关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。计算方法是有区别的。第二，误差中的均方差第二，误差中的均方差是未知的，必须使用其估计值是未知的，必须使用其估计值2112)1(1NiiNiiXNXN来代替，在计算所求量的同时，可计算出来代替，在计算所求量的同时，可计算出 。 7减小方差的各种技巧 显然，当给定置信度显然，当给定置信度后，误差后，误差由由和和N N决定。决定。
6、要减小要减小，或者是增大，或者是增大N N，或者是减小方差，或者是减小方差2 2。在。在固定的情况下，要把精度提高一个数量级，试验次固定的情况下，要把精度提高一个数量级，试验次数数N N需增加两个数量级。因此，单纯增大需增加两个数量级。因此，单纯增大N N不是一个不是一个有效的办法。降低方差的各种技巧，引起了人们的有效的办法。降低方差的各种技巧，引起了人们的普遍注意。普遍注意。一般来说，降低方差的技巧，往往会使观察一个子样一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观
7、察一个子样所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中特卡罗方法中效率效率的概念。的概念。它定义为 其中其中c c是观察一个子样的平均费用。是观察一个子样的平均费用。c283.蒙特卡罗方法的特点蒙特卡罗方法的特点 优点优点1)1)能够比较逼真地描述具有能够比较逼真地描述具有随机性质的事物的特点及随机性质的事物的特点及物理实验过程。物理实验过程。2)2)受几何条件限制小。受几何条件限制小。3)3)收敛速度与问题的维数无收敛速度与问题的维数无关。关。4)4)误差容易确定。误差容易确定。5)5)
8、程序结构简单，易于实现。程序结构简单，易于实现。 缺点缺点1)收敛速度慢。收敛速度慢。2)误差具有概率性。误差具有概率性。9第二章 随机数的产生1( )inf :( ),01Fyx F xyy2.1 2.1 逆变换法逆变换法设随机变量设随机变量X X的分布函数为的分布函数为F(x),F(x),定义定义定理定理2.1 设随机变量设随机变量U服从服从U(0,1)分布，则分布，则的分布函数为的分布函数为F(x).由定理由定理2.1,要生成分布函数为要生成分布函数为F(x)的随机数，可先生的随机数，可先生成成U(0,1)随机数随机数U，则可得到随机数，则可得到随机数X=F-1(U) 1XFU102.2
9、 合成法合成法如果如果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的密度函数为的密度函数为 110,1nniiiiiip xpx其中，由合成法，由合成法，X的随机数可如下抽取：的随机数可如下抽取：1）取）取uU(0,1); 2)取取 ,确定确定i，使，使3) 由由pi(x)抽取抽取
10、x.00100iijjjju112.3 筛选抽样筛选抽样 当当p(x)难以直接抽样时，如果可以将难以直接抽样时，如果可以将p(x) 表示成表示成p(x)=c h(x)g(x),其中其中h(.)是一密度函数且易于抽样，是一密度函数且易于抽样，而而0g(y),回到回到1）上述方法就是筛选抽样法，它是一种非常重要的抽样上述方法就是筛选抽样法，它是一种非常重要的抽样方法，可解决许多难以直接抽样的分布的抽样问题。方法，可解决许多难以直接抽样的分布的抽样问题。12h(x)的的选取有多种方法。一种直观的方法是：如果的的选取有多种方法。一种直观的方法是：如果存在一个函数存在一个函数M(x)，满足，满足p(x)
11、 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）。）。( )cM x dx 筛选抽样的理论依据如下：筛选抽样的理论依据如下：定理定理 设设X的密度函数为的密度函数为p(x),且且p(x)=c h(x)g(x)，其中，其中01时，如果时，如果 ，则，则x=y， 否则转到否则转到1）；）；yue1uy152.4 随机向量的抽样法随机向量的抽样法设设X1,Xk的联合概率密度为
12、的联合概率密度为11122111( ,.,)|.|,.,kkkkp xxpxpxxpxxx定理定理2.4 设设U1,Uk是独立同分布的是独立同分布的U(0,1)变量，变量， X1,Xk是方程是方程1111()|,.,2,.,iiiiF XUF XXXUik的解，其中的解，其中 是对应于是对应于 的分布函数，则的分布函数，则X1,Xk的的分布为分布为(2.4).iFip(2.4)(2.5)16随机向量的逆变换抽样法：随机向量的逆变换抽样法：1)由由U(0,1)分布独立地抽取分布独立地抽取u1,uk;2)用方程（用方程（2.5)解解x1,xk例例2.3 设设X1,X2的联合密度函数为的联合密度函数
13、为112121261,0,0( ,)0 xxxxxp x xothers试生成试生成X1,X2的随机数。的随机数。解：解：1111121110()66101xp xx dxxxx1222121111,1(|), 011p x xpxxxxpxx 17相应的边际分布函数和条件分布函数分别为相应的边际分布函数和条件分布函数分别为1231121110()6 132, 01xF xtt dxxxx21221212101(|)1, 011xF xxdtxxxxt 方程方程(2.5)变为变为231111212321XXUXXU此方程不易解，不妨交换两此方程不易解，不妨交换两自变量的次序自变量的次序2122
14、211220()63 101xpxx dxxx212112121222,(|)21, 01p x xp xxxxxxpx 18相应的边际分布函数和条件分布函数分别为相应的边际分布函数和条件分布函数分别为12122222221011121201()111;(|)11xxxF xxxF xxxx 方程方程(2.5)变为变为1223112111212212232121111111XUXXUXXUXUU 对服从特定分布的随机向量有一些特殊的抽样方法。对服从特定分布的随机向量有一些特殊的抽样方法。19例例2.6 试生成试生成k维正态分布维正态分布 的随机数。的随机数。,N解：注意到若解：注意到若 ，则存
15、在下三角阵，则存在下三角阵0,kZNI112122111111000cccCccc使使,kXCZN其中其中C可由迭代实现：可由迭代实现：首先，由首先，由 ，有，有11111Xc Z211111()Var Xc从而从而1111c。因。因22112112Xc Zc Z于是于是22222212212121221(),(,)Var XccCov XXc c得得212122122221111,.cc依此类推，依此类推，20一般迭代公式为一般迭代公式为11121,1,., ,1,.,jijiljllijjjjjllc ccik jic至此，我们可以给出至此，我们可以给出k维正态分布的抽样步骤：维正态分布的
16、抽样步骤：1)迭代计算迭代计算 ；2)由由N(0,1)分布独立抽取分布独立抽取k个随机数个随机数 ；3)计算计算,1,., ,1,.,ijc ik ji1,kzzzxCz212.5 随机模拟计算随机模拟计算2.5.1 随机投点法随机投点法 baf x dx考虑积分考虑积分 ，设，设a,b有限，有限，0 f(x) M,令令=(x,y):a=(x,y):a x x b,0b,0 y y M,M,并设并设(X,Y)(X,Y)是在是在上均匀上均匀分布的二维随机向量，其联合密度函数为分布的二维随机向量，其联合密度函数为则易见，则易见， 是是中曲线中曲线f(x)f(x)下方面积下方面积。 假设我们向假设我
17、们向中投点，若点落在中投点，若点落在y=f(x)下方称为中的，下方称为中的，则点中的概率为则点中的概率为,01,1a x by Mp x yM ba baf x dx 1bapf x dxM ba22若我们进行了若我们进行了n次投点，其中次投点，其中n0次中的，则可以得次中的，则可以得到到一个估计一个估计 01baM ba nf x dxn不难看出，不难看出， 是是的无偏估计，且其方差为的无偏估计，且其方差为1 11VarM baO nn(2.5.1)232.5.2 样本均值法样本均值法于是，积分于是，积分 baf x dxba E f x注意到，若注意到，若XU(a,b),则则 baf xE
18、 fXdxba由大数定律，若由大数定律，若 ，则，则21nPiibafXn . .1,.,. ( , )i i dnXXU a bMC方法为：方法为：1) 独立产生独立产生n个个U(a,b)随机数随机数2)按按(2.5.2)估计估计。1,.,nXX(2.5.2)24可证，在可证，在0 f(x) M条件下，条件下， 212211niibabaVarVarfXnbafx dxVarn2.5.3 降低方差的技术降低方差的技术Monte Carlo 方法中一类重要的研究课题是考虑一方法中一类重要的研究课题是考虑一些降低估计方差的技术。常用的方法有：重要抽样些降低估计方差的技术。常用的方法有：重要抽样法
19、，分层抽样法，关联抽样法等。法，分层抽样法，关联抽样法等。一一 重要抽样法重要抽样法由上节，样本平均法比投点法有效，将样本平均法做由上节，样本平均法比投点法有效，将样本平均法做更一般的推广，设更一般的推广，设g(x)是是(a,b)上的密度函数，改写上的密度函数，改写25由大数定律，若由大数定律，若 ，则，则311niPiifXfXEng Xg X . .1,.,. ( )i i dnXXg xMC方法为：方法为：1)选择适当的选择适当的g(x),独立产生独立产生n个个g(x)随机数随机数2)由由(2.5.3)估计估计。显然显然1,.,nxx 32minminggfXVarVarVarg X b
20、af xfXg x dxEg xg X(2.5.3)26从理论上看，因 2221fXVarEng X，若f(x)0,取 f xg x则有 20Var因为未知，这是作不到的，但它提示我们取g(x)与f(x)形状接近，应能降低方差。这就是重要抽样法的基本思想。其方差与其方差与g(x)有关。问题变为，如何选择有关。问题变为，如何选择g(.)使估使估计的方差最小。计的方差最小。27例例2.5.1 分别用投点法，均值法，重要抽样法，求积分别用投点法，均值法，重要抽样法，求积分分 ,比较各种方法的有效性。比较各种方法的有效性。10 xe dx解解 i)投点法投点法1)产生随机数 2) 对每对 ，记 的次数
21、为n0.则G(0,1),0,1,., ;iiXUYUein,iiX YiXiYe 01111.718,eenVarennnnii)均值法均值法1)产生随机数 2)(0,1),1,., ;iXUin 1222221011130.242222inXxieeVare dxennnn28iii)重要抽样法重要抽样法由重要抽样法的思想，需选择一个与 相似的密度函数。由Taylor展开式 取xe212xxex 1)产生随机数 2) 取则012( )113xg xx (0,1),1,., ;iuUin 131 1iiiXFuu 312122230132 113130.026912 12 1iXniiXxen
22、XeeVarEdxenXnxn（数值计算）真值真值投点法投点法均值法均值法重要抽样法重要抽样法1.718281.87561.81451.7219模拟结果29二、分层抽样法二、分层抽样法另一种利用贡献率大小来降低估计方差的方法是分层另一种利用贡献率大小来降低估计方差的方法是分层抽样法。它首先把样本空间抽样法。它首先把样本空间D分成一些不交的小区间分成一些不交的小区间 ，然后在各小区间内的抽样数由其贡献大小，然后在各小区间内的抽样数由其贡献大小决定。即，定义决定。即，定义 ，则，则Di内的抽样数内的抽样数ni应与应与pi成正比。成正比。1miiDD iiDpfx dx考虑积分考虑积分 10f x
23、dx将将0,1分成m个小区间：010.1maaa则则 11110iiammmiiaf x dxf x dxI记记 为第为第i个小区间的长度，个小区间的长度，i=1,m.在每个在每个小区间上的积分值可用均值法估计出来，然后将其小区间上的积分值可用均值法估计出来，然后将其相加即可给出相加即可给出的一个估计。具体步骤为：的一个估计。具体步骤为：1iiilaa301) 独立产生独立产生U(0,1)随机数随机数2)计算计算3) 计算计算,1,1,ijiujn im1,1,1,ijiiijixalujn im1iliiijjilIf xn于是可得于是可得的估计为的估计为41miiI(2.5.4)易见，易见
24、， 是是的无偏估计，其方差为的无偏估计，其方差为4 1224111222,iiinmmiiiijijiiaiiiiallVarVarfXnnfxIdxll 其中(2.5.5)(2.5.6)31续例续例2.5.1 考察分层抽样法求积分考察分层抽样法求积分 的方差。的方差。解：先将区间解：先将区间0,1划分成两个小区间划分成两个小区间0,0.5,0.5,1,则则10 xe dx0.511200.50.522210122220.51;2410.03492;240.09493;xxxxIe dxeIe dxeee dxee dxee设一共抽设一共抽n个随机数，其中在个随机数，其中在0,0.5)上抽上抽
25、n1个，则个，则使用分层抽样法求得使用分层抽样法求得 的方差为的方差为4 2222412110.50.5Varnnn32对对n1求导易知，在求导易知，在n固定下，当固定下，当 时时11120.37753nn4的方差最小，为的方差最小，为 420.06125VarVarn如果我们将区间进行如果我们将区间进行10等份，并确定出最优的抽样等份，并确定出最优的抽样次数分配：次数分配： ，则可得到分层抽样法，则可得到分层抽样法估计的方差为估计的方差为.1miijjn n 30.00246 nVar一般地，若诸一般地，若诸 已知，在已知，在n固定下，当固定下，当时，估计的方差最小，为时，估计的方差最小，为
26、,iil1miiiiiinnll211miiiln33分层抽样法在实施上有两个主要问题，其一是怎样划分区间，简单而常用的方法是将区间等分；另一个问题是在区间划分好后如何确定抽样次数的分配。由于在实际中 总是未知的，因而前面最优分配的结论无法应用。即使如此，分层抽样法还是有其作用的。可以证明，即使取简单的分配也有iiinnlba 42VarVar事实上，取 ，代入(2.5.5)得iinnlba 241miiibaVarln由Cauchy-Schwarz不等式，有2222211111mmmmmiiiiiiiiiiiiiiIIIIllballl据此，在(2.5.6)式两端各乘以 并相加得 il 22
27、222211bbmmiiiiiiaaIlfx dxfx dxlba于是于是 42VarVar34三、关联抽样法三、关联抽样法考虑积分差考虑积分差 1212fx dxfx dxII若用若用 估计估计，则其方差为，则其方差为12II 121212 2,VarVar IIVar IVar ICov I I显然，在显然，在 确定后，确定后， 正相关度越高，正相关度越高，则则 的方差越小。这便是关联抽样法的基本出发点。的方差越小。这便是关联抽样法的基本出发点。 11Var IVar I和12 ,I I考虑用重要抽样法来估计考虑用重要抽样法来估计I1,I2,即改写即改写为为 1122h x gx dxhx
28、 gx dx1)产生产生n个个U(0,1)随机数随机数 ；令；令2)则则1,nuu1112,iiiixGuyGu 1211niiih xhyn35第三章 数据添加算法 在Bayes统计或极大似然估计的计算中，经常会遇到这样一类问题：设我们能观测到的数据是Y，关于Y的后验分布p(|Y)很复杂，难以直接进行各种统计计算.假如我们能假定一些没有能观察到的潜在数据Z为已知（譬如，Y为某变量的截尾观测值，Z为该变量的真值），则可能得到一个关于的简单的添加后验分布p(|Y,Z),利用p(|Y,Z)的简单性我们可以进行各种计算，如极大化，抽样等，然后回过头来，又可以对Z的假定做检查或改进。如此进行，我们就将
29、一个复杂的极大化问题转变为一系列简单的极大化或抽样。在统计上，这种处理问题的方法称为“数据添加算法”。 常用的“数据添加算法”有EM算法和Markov Chain Monte Carlo方法。363.1 EM算法 先考虑一种简单情形。设某元件的失效时间Y关于变量x有直线回归关系，假设在一次试验中得到一批数据，如图， “ ”表示该元件失效时间坐标， ”“表示对应元件的截尾时间（小于失效时间）。如果直线斜率和截矩的估计值已知，则我们可以在真实数据不小于截尾数据的前提下将各个被截尾的失效时间估计出来，从而得到所谓的”完全数据“，由此完全数据，重新对直线的斜率及截矩进行估计，再依据新的估计量，得到新的
30、”完全数据“。如此循环往复，则将一个复杂的估计问题替换成一系列简单的估计问题。将之一般化，就给出EM算法。37EM算法是一种迭代方法，主要用来求后验分布的众算法是一种迭代方法，主要用来求后验分布的众数（即极大似然估计）。它的每一步迭代由两步组数（即极大似然估计）。它的每一步迭代由两步组成：成：E步（求期望）和步（求期望）和M步（极大化）。步（极大化）。一般地，以一般地，以p(|Y)|Y)表示表示基于基于Y Y的的后验密度，称为的的后验密度，称为观测后验分布；观测后验分布； p(|Y,Z)表示添加数据表示添加数据Z后得到的后得到的的后验密度，称为添加后验分布；的后验密度，称为添加后验分布； p(
31、Z|,Y)表示在给表示在给定观测数据定观测数据Y和参数和参数条件下条件下Z Z的条件密度。的条件密度。我们的目的是计算我们的目的是计算p(|Y)的众数。的众数。于是于是EM算法如下进行。记算法如下进行。记 为第为第i+1次迭代开始时次迭代开始时后验众数的估计值，则第后验众数的估计值，则第i+1次迭代的两步为次迭代的两步为E步：将步：将p(|Y,Z)或或log p(|Y,Z)关于关于Z的条件分布求期的条件分布求期望，从而把望，从而把Z积掉，即积掉，即 i38 |,log| ,|,log| ,|,iiizQYEpY ZYpY z p zY dZ M步：将步：将 极大化，即找到一个点极大化，即找到一
32、个点 ，使，使 |,iQY 1i 1|,max|,iiiQYQY 将上述将上述E，M步循环进行，直至步循环进行，直至 充分小为充分小为止。止。 1ii例例3.1 设总体设总体X的分布律为的分布律为其中其中 (0,1),(0,1),现进行了现进行了X1234pk1241141144197次试验，观察到次试验，观察到1,2,3,4的频数为的频数为取取的先验分布的先验分布( () )为为U(0,1)U(0,1)分布分布, ,则则的观察后的观察后验分布为验分布为125,18,20,34Y 39 1234123411|1244421yyyyyyyypYp Y 现假设现假设X=1可以分解为两部分，其发生概
33、率分别为可以分解为两部分，其发生概率分别为1/2和和/4,/4,令和令和y y1 1-Z-Z分别表示试验结果中落入这分别表示试验结果中落入这两部分的次数（是不能观测到的潜在数据），则两部分的次数（是不能观测到的潜在数据），则的添加后验分布为的添加后验分布为 1234231411|,|124441zyzyyyyyyzypY Zp Y Z (3.1.1)(3.1.2)显然，用显然，用(3.1.2)式求极值比式求极值比(3.1.1)式简单。迭代如下：式简单。迭代如下：40E步：在给定下，步：在给定下， 1423142311423|,loglog 1|,|,loglog 12loglog 12iizi
34、ziQYEyzyyyYyEzYyyyyyyyy ,iY 12,2iZb yM步：将步：将 关于关于极大化得极大化得 |,iQY 14112342159681971442iiiiiiiyyyyyy可以证明，在关于logp(|Y)的很一般的条件下，由算法得到的估计序列 收敛到的稳定点。（不能保证是极大值点）。较为可行的办法是选几个不同的初值迭代，然后在诸估计值中加以选择，这可减轻初值选取对结果的影响） i41估计的精度假设假设EM算法最后的结果是算法最后的结果是 ，则根据似然估计的渐，则根据似然估计的渐近正态性，其渐近方差可用近正态性，其渐近方差可用 Fisher观测信息的倒数观测信息的倒数121
35、002log|pYVI 近似。近似。（证明见高等数理统计p126定理2.5.4)423.2 Markov Chain Monte Carlo方法 ,xxX Eff xx dxX对于较简单的后验分布，可直接计算或静态对于较简单的后验分布，可直接计算或静态MC等近等近似计算方法。但在实际中，观测后验分布往往是复杂似计算方法。但在实际中，观测后验分布往往是复杂的，高维的，非标准形式的分布，上述方法都难以实的，高维的，非标准形式的分布，上述方法都难以实施。对于这类问题，一种简单且行之有效的施。对于这类问题，一种简单且行之有效的Bayes计计算方法就是算方法就是MCMC。EM算法得到的是后验分布的众数，
36、有时我们希望算法得到的是后验分布的众数，有时我们希望得到其它一些后验量如后验均值，方差，后验分布得到其它一些后验量如后验均值，方差，后验分布的分位数等。计算这些后验量都可归结为关于后验的分位数等。计算这些后验量都可归结为关于后验分布积分的计算。具体地，设分布积分的计算。具体地，设 为后验密为后验密度，我们要计算的后验量可写成某函数度，我们要计算的后验量可写成某函数f(x)关于关于的的期望期望(3.2.1)433.2.1 基本思路MCMC方法的基本思想是通过建立一个平稳分布为方法的基本思想是通过建立一个平稳分布为(x)(x)的的MarkovMarkov链来得到链来得到(x)(x)的样本，基于这些
37、样本的样本，基于这些样本可以作各种统计推断。比如，若得到了平稳分布为可以作各种统计推断。比如，若得到了平稳分布为(x)(x)的的Markov链的链的样本轨道样本轨道 ，则，则(3.2.1)(3.2.1)可估计为可估计为 1,.,nXX 1nimni mffXnm(3.2.2)注注 由Markov链平稳分布的概念可知，不论Markov链从什么初始状态出发，经过一段时间后，各个时间的边际分布都是平稳分布，因此可将经过某个m时间之后的观察值看作平稳分布(x)的样本。由遍历性定理可知，,mnfE f n MCMC的关键是如何构造平稳分布为的关键是如何构造平稳分布为的的Markov链的转移核p(x,y)
38、44MCMC方法可概括为如下三步：方法可概括为如下三步：(1) 在在X X上选一个上选一个“合适合适”的的Markov链，确定其转链，确定其转移核移核p(x,y),使链的平稳分布为使链的平稳分布为。(2)(2)由由X X中某一点中某一点X X(0)(0)出发，用出发，用(1)(1)中的中的Markov链产链产生序列生序列X1,Xn;(3) 对某个对某个m和大的和大的n，任一函数，任一函数f(x)的期望估计如下的期望估计如下11nii mE ffXnmMCMC有许多研究专题，如链的收敛性判断（m大小的确定），链的长度（n的大小）的确定，估计误差等等。以下主要讨论转移核的构造。453.2.2 满条
39、件分布MCMC主要用于多变量，非标准形式，且各变量间主要用于多变量，非标准形式，且各变量间相互不独立时分布的模拟。相互不独立时分布的模拟。令令 ，我们总可以写出，我们总可以写出1,.,nxxx 1|niiixxx其中 。如果(3.2.1)式中右端各个因子能够直接模拟，则只需要进行静态模拟（抽样过程中不改变抽样分布）。实际中很难满足上述条件，因此需进行动态模拟（抽样分布随模拟的进行而改变），如MCMC，此时满条件分布扮演了一个重要角色。,ijxxji(3.2.1)在导出满条件分布时，应注意到这样一个事实：记 ，x X,ijxxji,ijxxji,TjxxjT |TTTxxxxx dx(3.2.2
40、)46等价地，若等价地，若 ，且，且 ，则，则, x x XTTxx |TTTTxxxxxx(3.2.3)一般地，用一般地，用y表示观测数据，表示观测数据， ，其中，其中 分分别表示参数，超参数和缺损数据，则有别表示参数，超参数和缺损数据，则有, ,xz , ,z |, , , |x yy zp y z 其中，其中， 表示完全数据的密度函数，表示完全数据的密度函数， 表表示先验分布，示先验分布， 表示超参数的分布。有表示超参数的分布。有(3.2.2),各变各变量的满条件分布如下：量的满条件分布如下：( , | )p y z| |, , , |,iiiiz yp y z | , ,|,jjz y
41、 | , ,( , | )kkzzyp y z 47例例3.2.1 设设(X1,X2)的联合密度为的联合密度为2211 222122(1)1221( ,),21xx xxx xe且且 ，则其满条件分布为，则其满条件分布为0,1U22222122(1)01222(1)0121|,1111,1,1,2,2,11iijjijxx xxijxxjxxeeNxij 2211222122 (1)1221|,1xx xxxxe483.2.3 Gibbs 抽样思想：设思想：设 的密度为的密度为 ，任意固，任意固定定T N，在给定，在给定 条件下，如下定义随机变条件下，如下定义随机变量量 具有密度函具有密度函数
42、数 ，则对任一可测集，则对任一可测集B，1,nXXX xTTXx1 ,:,nTTTXXXXXX而 |TTxx |TTTBBP XBxxxdxx dxB因而因而X的密度也是的密度也是 。 x上述过程定义了一个由上述过程定义了一个由X到到X的转移核，且其相应的转移核，且其相应的平稳分布是的平稳分布是。这样构造的。这样构造的MCMCMCMC称为称为GibbsGibbs抽样。抽样。当当T T只有一个元素时称为单元素只有一个元素时称为单元素Gibbs抽样。抽样。49单元素单元素Gibbs抽样具体步骤如下：抽样具体步骤如下：在给定起始点在给定起始点 后，假定第后，假定第t次迭代开次迭代开始时的估计值为始时
43、的估计值为 ，则第，则第t次迭代分为如下次迭代分为如下n步：步：（1）由满条件分布）由满条件分布 抽取抽取 (i)由满条件分布由满条件分布 抽取抽取 ; (n)由满条件分布由满条件分布 抽取抽取 0001,.,nxxx1tx1112|,.,ttnxxx 1tx 11111|,.,.,ttttiiinxxxxx tix 11|,.,ttnnxxx tnx记记 ，则，则 是平稳分布为是平稳分布为的的Markov链的实现值，其由链的实现值，其由x到到x的转移概率函数为的转移概率函数为 1,.,tttnxxx 12,txxx 1221311, |,., | ,.,. | ,., nnnnp x xxx
44、xxxxxxxx503.2.3 Metroplis-Hastings方法|iixx在在Gibbs抽样中，抽样中， 可能很难抽取，这时可采用可能很难抽取，这时可采用更一般化的更一般化的Metroplis-Hastings方法。其思路如下：方法。其思路如下： 任意选择一个不可约的转移概率任意选择一个不可约的转移概率q(.,.)以及一个以及一个函数函数(.,.)(.,.)，00 1,1,对任一组合对任一组合(x,x(x,x)(x)(x x x),),定定义义 , , , p x xq x xx xxx则则p(x,x)形成一个转移核。形成一个转移核。此方法的实施比较直观，首先由此方法的实施比较直观，首
45、先由q(x,x)产生一个潜在产生一个潜在的转移的转移x x，然后根据概率，然后根据概率(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。51注意：在有了注意：在有了q(.,.)q(.,.)后，应选择后，应选择(.,.)(.,.)使相应的使相应的p(x,xp(x,x) )以以(x)为平稳分布。为平稳分布。定理定理3.2.1 3.2.1 给定给定q(x,xq(x,x),),取取 , min 1,
46、 xq x xx xx q x x则由以下则由以下p(x,x)产生的产生的Markov链是可逆的，且以链是可逆的，且以(x)为平稳分布为平稳分布 , , , , , , , p x xq x xx xq x xx q x xx q x xxq x xx q x xx q x xx52由定理由定理3.2.1可以看出，建议分布可以看出，建议分布q(x,x)可取各种形式。可取各种形式。常用的建议分布选择方法：常用的建议分布选择方法：(1)Metropolis选择：对称的建议分布，即选择：对称的建议分布，即,q x xq x x此时此时 , min 1,xx xx例如例如 22exp,exp22xxx
47、xq x x53（2）独立抽样）独立抽样取取 ,q x xq x则则 , min 1,xq xx xx q x一般，要使独立抽样有好的效果，一般，要使独立抽样有好的效果，q(x)应接近应接近(x),(x),比较安全的办法是使比较安全的办法是使q(x)q(x)的尾比的尾比(x)重。重。(3)单元素单元素Metropolis-Hastings算法算法产生向量产生向量X的样本有时是困难的，通常是利用满的样本有时是困难的，通常是利用满条件分布对条件分布对X的分量进行逐个抽样。的分量进行逐个抽样。考虑考虑 的条件分布，选择一个转移核的条件分布，选择一个转移核 ，固定，固定 不变，由不变，由q 产生一产生
48、一个可能的个可能的 ，然后以概率，然后以概率决定是否接受其作为下决定是否接受其作为下一状态。一状态。|,1,.,iiXXin|iiiq xxxiiiXXx ix54（续）例（续）例3.2.1 设设(X1,X2)的联合密度为的联合密度为2211 222122(1)1221( ,),21xx xxx xe试产生试产生 的后验分布样本。的后验分布样本。解：由例解：由例3.2.1, 各分量的满条件分布为各分量的满条件分布为2|,1,1,2,2,1ijjxxNxij2211222122 (1)1221|,1xx xxxxe55 2211 2222211 222122(1)2122(1)211,min 1,min 1,11xx xxxx xxee 1212|,|,1,(01)qx xqx x取取则抽样步骤为：则抽样步骤为：(1)选初值选初值(2) 按按 产生随机数产生随机数(3)按按 U(0,1)产生随机数产生随机数(4)(4)产生一产生一U(0,1)U(0,1)随机数随机数u,u,如果如果 则否则否则则 ，否则，否则 0010,1 ,xR 000 21,1Nx 02x 01,u 1 10此课件下载可自行编辑修改，供参考！此课件下载可自行编辑修改，供参考！感谢您的支持，我们努力做得更好！感谢您的支持，我们努力做得更好！
