白 第11章Markov链Monte Carlo方法 §11.1 计算积分的Monte Carlo方法 §11.2 Markov链Monte Carlo方法简介 §11.3 Metropolis-Hastings算法 §11.4 Gibbs抽样 §11.5 贝叶斯MCMC估计方法 GoBack FullScreen Close Quit
2/52 kJ Ik J I GoBack FullScreen Close Quit 111Ÿ MarkovÛMonte Carloê{ • § 11.1 O黩Monte Carloê{ • § 11.2 MarkovÛMonte Carloê{{0 • §11.3 Metropolis-Hastingsé{ • §11.4 Gibbsƒ • §11.5 ìdMCMCOê{
Monte Carlo方法也即随机模拟方法的别称，在应用 统计学、金融学、经济计量学、计算物理学等多个领域有着 广泛的应用。法国数学家蒲丰（Buffon)于1777年提出的 利用投针实验问题，利用针与平行线相交的频率近似表示 相交概率，从而求得圆周率，这被看成是Monte Carlo方 法的起源。而由“计算机之父”和“博弈论之父”冯.诺依 曼(John von Neumann)等制订的使用随机数处理确定 3/52 性数学问题的Monte Carlo方法极大推动了数值分析的发 展。冯.诺依曼使用摩洛哥的赌城Monte Carlo来命名这种 方法。 GoBack FullScreen Close Quit
3/52 kJ Ik J I GoBack FullScreen Close Quit Monte Carloê{è=ëÅ [ê{O°ß3A^ ⁄OÆ!7KÆ!²LO˛Æ!Oé‘nÆ ıá+çkX 2çA^"{IÍÆ[Æ¥£Buffon§u1777cJ— |^›¢ØKß|^ܲ1ÇÉ™«CqL´ ÉV«ßl ¶ ±«ß˘w§¥Monte Carloê { " d/OéÅÉI0⁄/ÆâÿÉI0æ.Ïù ˘£John von Neumann§ õæ¶^ëÅÍ?n(½ 5ÍÆØKMonte Carloê{4å̃ Í䩤u –"æ.Ïù˘¶^ ‚xŸ¢Monte Carlo5·¶˘´ ê{"
Monte Carlo方法的基本原理是：当求解随机事件发 花 生的概率或随机变量的数学期望时，通过设计的某种实验， 得出某个特定事件发生的频率，使用这个频率来近似表示 这一事件发生的概率，从而得到问题的数值解。在Monte Carlo方法中包含三个核心问题，分别为：构造概率过程、 从已知概率分布中抽样、建立估计量。 4/52 GoBack FullScreen Close Quit
4/52 kJ Ik J I GoBack FullScreen Close Quit Monte Carloê{ƒ n¥µ¶)ëÅØáu )V«½ëÅC˛ÍÆœ"ûßœLO,´¢ß —,áA½Øáu)™«ß¶^˘á™«5CqL´ ˘òØáu)V«ßl ØKÍä)"3Monte Carloê{•ù¹náÿ%ØKß©OèµEV«Lß! lÆV«©Ÿ•ƒ!Ô·O˛
y 通过构造独立同分布的随机数来计算积分的Monte Carlo方法，也称为静态Monte Carlo方法。在维数非常 高的情况下，因为计算量太大，使用静态Monte Carlo方 法处理速度太慢。故引入一种动态Monte Carlo,即通过 构造Markov链的极限平稳分布，来模拟计算积分的方法， 称为Markov链Monte Carlo方法，以下均简记为MCMC (Markov Chain Monte Carlo)。本章简要介绍计算积 5/52 分的Monte Carlo方法，以及Metropolis-Hastings.算法， Gibbs抽样方法。 GoBack FullScreen Close Quit
5/52 kJ Ik J I GoBack FullScreen Close Quit œLE’·”©ŸëÅÍ5O黩Monte Carloê{ßè°è·Monte Carloê{"3ëÍö~ pú¹eßœèOé˛å߶^·Monte Carloê {?nÑ›˙"⁄\ò´ƒMonte Carloß=œL EMarkovÛ4Ų©Ÿß5 [O黩ê{ß °èMarkovÛMonte Carloê{ß±e˛{PèMCMC £Markov Chain Monte Carlo§" Ÿ{á0 Oé» © Monte Carloê{ß±9Metropolis-Hastingsé{ß Gibbsƒê{
§11.1 计算积分的Monte Carlo方法 静态Monte Carlo方法是通过构造独立同分布的随机 数来计算积分，有频率法与期望法两种，我们分别给出例 子。其中频率法以经典的Buffon投针问题为例。 例11.1.1平面上画着一些平行线，线之间的距离均 为a，向此平面随机地投一长度为(l<a)的针投掷到平面 6/52 上，求针与任一平行线相交的概率。 解：以Z表示针的中点到最近一条平行线的距离，B表 示针与平行线的交角，则(B,Z)可以确定针所落的位置，且 有0≤Z≤fraca2,0≤B≤π，即针的所有可能位置可 用BZ平面上的矩形来表示，将此矩形区域记为G。 GoBack FullScreen Close Quit
6/52 kJ Ik J I GoBack FullScreen Close Quit §11.1 O黩Monte Carloê{ ·Monte Carloê{¥œLE’·”©ŸëÅ Í5O黩ßk™«{Üœ"{¸´ß·Ç©Oâ—~ f"Ÿ•™«{±²;Buffon›ØKè~" ~ 11.1.1 ²°˛xXò ²1ÇßÇÉmÂl˛ èa ßïd²°ëÅ/›ò›èl(l < a)›ï²° ˛ß¶Ü?ò²1ÇÉV«" )µ±ZL´•:ÅCò^²1ÇÂlßβL ´Ü²1ÇßK(β, Z)å±(½§·†òßÖ k0 ≤ Z ≤ fraca2, 0 ≤ β ≤ πß=§kåU†òå ^βZ²°˛›/5L´ßÚd›/´çPèG"
父 当针与平行线相交时，Z与B应满足关系Z≤号sm,B,即 为图11-1(b)中的阴影部分，将此区域记为R. a/2 [sinB/2 G R π B 7/52 图11-1 Buffon投针图 故针与平行线相交的概率为： L(R) IsinBdB 21 p= L(G) ira πa GoBack FullScreen Close Quit
7/52 kJ Ik J I GoBack FullScreen Close Quit ܲ1ÇÉûßZÜβA˜v'XZ ≤ l 2 sinβß= è„11-1(b)•“K‹©ßÚd´çPèR. „11-1 Buffon›„ ܲ1ÇÉV«èµ p = L(R) L(G) = 1 2 R π 0 lsinβdβ 1 2 πa = 2l πa
利用此例的答案来估算π的值，其方法是投针N次，记 录针与线相交的次数，再以频率值朵作为概率p的值，于 是可得 2l N ≈ πa 即得 2IN π≈ an 8/52 以下我们将使用R软件来实现Buf fon投针实验，其中两个 关键的步骤为： (1)产生随机数：Z在[0，引内随机取值，B在[0，π]内随 机取值，也即首先需要生成服从[0，引上均匀分布的随机数， 以及服从0，π]上均匀分布的随机数。 GoBack FullScreen Close Quit
8/52 kJ Ik J I GoBack FullScreen Close Quit |^d~âY5éπäߟê{¥›NgßP ¹ÜÇÉgÍnß2±™«än NäèV«p äßu ¥å n N ≈ 2l πa = π ≈ 2lN an ±e·ÇÚ¶^R^á5¢yBuf f on›¢ßŸ•¸á 'Ö⁄½èµ £1§)ëÅ͵Z3[0, a 2 ]SëÅäßβ3[0, π]Së Åäßè=ƒkIá)§—l[0, a 2 ] ˛˛!©ŸëÅÍß ±9—l[0, π]˛˛!©ŸëÅÍ"
女 (2)模拟投针实验：考察次投针试验，若针与平行 线相交，则表示试验成功，假设共成功次，则π的估计值 为心，其中线之间的距离a及针的长度l(l<a)均预先给 an 定。 使用R软件来进行模拟，例如选取试验次数为N 50000,线之间的距离a=1,针的长度1=0.75，代码如 下： 9/52 GoBack FullScreen Close Quit
9/52 kJ Ik J I GoBack FullScreen Close Quit £2§ [›¢µ ng›£ßeܲ1 ÇÉßKL´£§ıßb§ıngßKπOä è2lN an ߟ•ÇÉmÂla9›l(l < a)˛˝kâ ½" ¶^R^á5?1 [ß~X¿£gÍèN = 50000ßÇÉmÂla = 1ß›l = 0.75ßìËX eµ
液 touzhen <-function(N,I=0.75,a =1) +{ +n<-0;beta <-runif(N,0,pi);z<-runif(N,O,a/2) +for(iin1:N) +{ 10/52 +if(z[i]<=l sin(beta[i])/2) +m<-n+1 +} +2*l*N/(n*a) +} touzhen(50000) GoBack FullScreen Close Quit
10/52 kJ Ik J I GoBack FullScreen Close Quit > touzhen touzhen(50000)
。 例11.1.2 (随机投点法)计算任一定义域为0,1]的 函数g(c)在区间[0,1]上的积分g(x)dxc. 解假设随机变量X和Y均服从[O,1上的均匀分布，且相 互独立，则二维随机向量(X,Y)服从矩形区域[0,1]*[0,1]上 的二维均匀分布，联合概率密度为： f(x,y）= 1,0<x<1,0<y<1 0, 11/52 其他 现用B表示事件{w:Y≤g(X)}(g为任-定义域为[0,1oe), 也即我们向矩形区域[0,1]*[0,1随机投点，其中点落在 以[0,1为底，以函数g(x)为曲边的曲边梯形内。 可计算事件B发生的概率为 P=PY≤9(X}-a f(x,y)dxdy GoBack FullScreen Close Quit
11/52 kJ Ik J I GoBack FullScreen Close Quit ~ 11.1.2 £ëÅ›:{§Oé?ò½¬çè[0, 1] ºÍg(x)3´m[0, 1]˛»© R 1 0 g(x)dx. ) bëÅC˛X⁄Y ˛—l[0, 1]˛˛!©ŸßÖÉ p’·ßKëëÅï˛(X, Y )—l›/´ç [0, 1]∗[0, 1]˛ ë˛!©ŸßÈ‹V«ó›èµ f(x, y) = ( 1, 0 >F/S" åOéØáBu)V«è P(B) = P{Y ≤ g(X)} = Z Y ≤g(X) f(x, y)dxdy