Lecture 8:Markov Chain Monte Carlo Methods (一） (马尔科夫蒙特卡罗方法) 张伟平 Monday 26th October,2009
Lecture 8: Markov Chain Monte Carlo Methods (ò) (ÍâÅÑAk¤ê{) ‹ï² Monday 26th October, 2009
Contents 1 Markov Chain Monte Carlo Methods y 1.1 Introduction..........·.。....······· 1.1.1 Integration problems in Bayesian inference...... 2 1.l.2 Markov Chain Monte Carlo Integration·.·..·· ￥ 1.1.3 arkov Chain.·.··.··.·········· b l.2 The Metropolis-Hastings Algorithm.·...·-······· 13 1.2.1 Metropolis-Hastings Sampler 14 1.2.2 The Metropolis Sampler................22 1.2.3 Random Walk Metropolis 23 1.2.4 The Independence Sampler 33 1.3 Single-component Metropolis Hastings Algorithms 37 l.4 Application:Logistic regression············· ...39 Previous Next First Last Back Forward 1
Contents 1 Markov Chain Monte Carlo Methods 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Integration problems in Bayesian inference . . . . . . 2 1.1.2 Markov Chain Monte Carlo Integration . . . . . . . 4 1.1.3 Markov Chain . . . . . . . . . . . . . . . . . . . . . 8 1.2 The Metropolis-Hastings Algorithm . . . . . . . . . . . . . . 13 1.2.1 Metropolis-Hastings Sampler . . . . . . . . . . . . . 14 1.2.2 The Metropolis Sampler . . . . . . . . . . . . . . . . 22 1.2.3 Random Walk Metropolis . . . . . . . . . . . . . . . 23 1.2.4 The Independence Sampler . . . . . . . . . . . . . . 33 1.3 Single-component Metropolis Hastings Algorithms . . . . . 37 1.4 Application: Logistic regression . . . . . . . . . . . . . . . . 39 Previous Next First Last Back Forward 1
Chapter 1 Markoy Chain Monte Carlo Methods 1.1 Introduction MCMC(Markov Chain Monte Carlo)方法的一般理论框架可以参看Metropo- 1 lis et al..(1953)以及Hastings(1970),以及其他各种介绍MCMC的专著.本 节我们介绍这种方法的基本思想和应用， 注意在前面介绍Monte Carlo方法估计积分 g(t)dt, 是把此积分表示成一个对某个概率密度f(t)下的期望.从而积分估计问题转化 成从目标概率密度(t)中产生随机样本。 Previous Next First Last Back Forward
Chapter 1 Markov Chain Monte Carlo Methods 1.1 Introduction MCMC(Markov Chain Monte Carlo) ê{òÑnÿµeå±Îw Metropolis et al. (1953)±9 Hastings (1970), ±9Ÿ¶à´0 MCMC;Õ. !·Ç0 ˘´ê{ƒ gé⁄A^. 5ø3c°0 Monte Carloê{O»© Z A g(t)dt, ¥rd»©L´§òáÈ,áV«ó›f(t)eœ". l »©OØK=z §l8IV«ó›f(t)•) ëÅ . Previous Next First Last Back Forward 1
在MCMC方法中，首先建立一个马尔科夫链，使得f(t)为其平稳分布.则可 以运行此马尔科夫链充分长时间直至收敛到平稳分布，那么从目标分布f(t)中 产生随机样本，就是从达到平稳状态的马尔科夫链中产生其样本路径.我 们将介绍几种建立这样的马尔科夫链的方法：Metropolis算法，Metropolis- Hastings算法，以及Gibbs抽样方法.一个好的链应该具有快速混合(rapid mixing)性质一从任意位置出发很快达到平稳分布. 1.1.1 Integration problems in Bayesian inference Bayesian推断中的许多问题都是MCMC方法的应用.从Bayesian的观点来看， 模型中的观测变量和参数都是随机变量.因此，样本x=(1,·,xm)和参 数联合分布可以表示为 fx,e(r,)=fxle(x1,…,xn)π(0). Previous Next First Last Back Forward 2
3MCMCê{•, ƒkÔ·òáÍâÅÛ, ¶f(t)蟲©Ÿ. Kå ±$1dÍâÅÛø©ûmÜñ ¬Ò²©Ÿ, @ol8I©Ÿf(t)• )ëÅ , “¥là²GÍâÅÛ•)Ÿ ¥ª. · ÇÚ0 A´Ô·˘ÍâÅÛê{: Metropolisé{, MetropolisHastingsé{, ±9Gibbsƒê{. òá–ÛAT‰kØÑ·‹(rapid mixing)5ü—l?ø†ò—uÈØ಩Ÿ. 1.1.1 Integration problems in Bayesian inference Bayesian̉•NıØK—¥MCMCê{A^. lBayesian*:5w, .•*ˇC˛⁄ÎÍ—¥ëÅC˛. œd, x = (x1, · · · , xn)⁄Î ÍθÈ‹©Ÿå±L´è fx,θ(x, θ) = fx|θ(x1, · · · , xn)π(θ). Previous Next First Last Back Forward 2
从而根据Byes定理，可以通过样本x=(x1,·,xm)的信息对0的分布进行更 新： f(x)x(0) fol(0x)== ∫fr1e(x)π(0)d0 则在后验分布下，9()的期望为 Eg()=/g(0)folz(0)de= ∫g(0)fzle(x)r(e)dB ∫fxla(x)π(e)d0 此积分为值为样本x的函数.因此可以对g(0)进行推断.比如g()=时， 则Eg(z)=E可以作为的估计. 对此类问题我们考虑一般形式： Eg(Y)= Jg(t)(t)dt ∫r(t)dt 这里π()为一个密度或者似然.若为一密度，则此期望即为常见的期望定义： Eg(Y)=∫9()fy(t)d.若π为一似然，则需要一个正则化常数才可以成为密 度.在贝叶斯分析中，π()为一后验密度.当π()的正则化常数未知时，此积分 Previous Next First Last Back Forward 3
l ä‚Bayes½n, 屜L x = (x1, · · · , xn)&EÈθ©Ÿ?1ç #: fθ|x(θ|x) == fx|θ(x)π(θ) R fx|θ(x)π(θ)dθ . K3©Ÿe, g(θ)œ"è Eg(θ|x) = Z g(θ)fθ|x(θ|x)dθ = R g(θ)fx|θ(x)π(θ)dθ R fx|θ(x)π(θ)dθ . d»©èäè xºÍ. œdå±Èg(θ)?1̉. 'Xg(θ) = θû, KEg(θ|x) = E[θ|x] å±äèθO. ÈdaØK·ÇƒòÑ/™: Eg(Y ) = R g(t)π(t)dt R π(t)dt . ˘pπ(·)èòá󛽈q,. eπèòó›, Kdœ"=è~Ñœ"½¬: Eg(Y ) = R g(t)fY (t)dt. eπèòq,, KIáòáKz~͂屧èó ›. 3ìd©¤•, π(·)èòó›. π(·) Kz~Íôû, d»© Previous Next First Last Back Forward 3
也可以计算.（分子分母中的正则化常数抵消了）.这是非常好的一个性质，因为 在很多问题下，正则化常数很难计算 另一方面，在很多具体问题中此积分根本没有显示表示，数值方法也很难计 算，特别是在高维时，而MCMC方法对此类积分提供了一个很好的计算方法. 1.1.2 Markov Chain Monte Carlo Integration 积分Eg(elx】=∫g()fa1r(r)d的Monte Carlo估计为样本均值 其中x1,·,工m为从分布faz(x)中抽取的样本.当x1,·,xm独立时，则可 以根据大数律知当样本量n趋向于无穷时，©收敛到Eg(x), 但是在一些问题中，从分布f1x(x)中抽取样本是非常困难的，MCMC方 法就是为此目的而诞生的，其第一个”MC",Markov Chain,就表示从目标分 布中抽样，第二个"MC",Monte Carlo,则表示在抽取到的样本下，应用Monte Previous Next First Last Back Forward 4
èå±Oé. (©f©1•Kz~Í-û ). ˘¥ö~–òá5ü, œè 3 ÈıØKe, Kz~ÍÈJOé. ,òê°, 3Èı‰NØK•d»©ä vkw´L´, Íäê{èÈJO é, AO¥3pëû, MCMCê{Èda»©J¯ òáÈ–Oéê{. 1.1.2 Markov Chain Monte Carlo Integration »© E[g(θ|x)] = R g(θ)fθ|x(θ|x)dθ Monte CarloOè ˛ä g¯ = 1 m Xm i=1 g(xi), Ÿ•x1, · · · , xmèl©Ÿfθ|x(θ|x)•ƒ . x1, · · · , xm’·û, Kå ±ä‚åÍÆ ˛n™ïuðû, ¯g¬ÒE[g(θ|x). ¥3ò ØK•, l©Ÿfθ|x(θ|x)•ƒ ¥ö~(J, MCMCê {“¥èd8 ), Ÿ1òá”MC”, Markov Chain, “L´l8I© Ÿ•ƒ, 1á”MC”, Monte Carlo, KL´3ƒ e, A^Monte Previous Next First Last Back Forward 4
Carlo积分方法对积分进行计算. MCM方法的理论依据在于下述的极限定理： Theorem1.设{X:,t≥0}为一不可约的，非周期的状态空间为2的马氏链， π为平稳分布，T0是初始分布，则 πt→r,t→00. 这里：表示马氏链在t时刻的边际分布.此定理说明了当此马氏链运行充分 长时间后(t=n,n很大)，则在时刻n,X的分布为霜，且和初始分布无关 Theorem2.(Markov chain Law of Large Numbers),若X1,X2,,·为一 历的平稳分布为π的马氏链值其中每个X:可以是多维的)，则X依分布收敛 到分布为π的随机变量X,且对任意函数f,当Ef(X)川存在时，且n→co时， 则 n=∑fK)→EfX).as m乡1 Previous Next First Last Back Forward 5
Carlo»©ê{È»©?1Oé. MCMê{nÿù‚3ue„4Žn: Theorem 1. {Xt, t ≥ 0}èòÿå,ö±œGòmèΩͺÛ, π販Ÿ, π0¥–©©Ÿ, K πt → π, t → ∞. ˘pπtL´ÍºÛ3tûè>S©Ÿ. d½n`² dͺÛ$1ø© ûm(t = n, nÈå), K3ûèn, Xn©Ÿèπ, Ö⁄–©©ŸÃ'. Theorem 2. (Markov chain Law of Large Numbers)eX1, X2, · · · èòH {²©ŸèπͺÛä(Ÿ•záXtå±¥ıë), KXnù©Ÿ¬Ò ©ŸèπëÅC˛X, ÖÈ?øºÍf, Eπ|f(X)|3û, Ön → ∞û, K f¯n = 1 n Xn t=1 f(Xt) → Eπf(X), a.s. Previous Next First Last Back Forward 5
下面我们计算遍历均值fn的方差.记f=f(X),Yk=cou(f,f+), 则∫的方差为a2=0,k步相关系数Pk=Y/a2.从而 店=awa=ava2刊=wr宫 =宫va:内+mr内 =+”=+ 可以证明 层→2=2+2∑p则 k=1 .因此遍历均值fn的方差为 n-1 k=1 Previous Next First Last Back Forward 6
e°·ÇOéH{˛äf¯nê . Pf t = f(Xt), γk = cov(f t , ft+k), Kf tê èσ 2 = γ0, k⁄É'XÍρk = γk/σ2 . l τ 2 n = nV arπ(f¯n) = nV ar( 1 n Xn t=1 f t ) = 1 n V arπ( Xn t=1 f t ) = 1 n Xn t=1 V arπ(f t ) + 2 1 n nX−t k=1 cov(f t , ft+k ) = σ 2 + 2 nX−1 k=1 n − k n γk = σ 2 [1 + 2 nX−1 k=1 n − k n ρk]. å±y² τ 2 n → τ 2 = σ 2 [1 + 2X∞ k=1 ρk] . œdH{˛äf¯nê è V arπ(f¯n) = σ 2 n [1 + 2 nX−1 k=1 n − k n ρk]. Previous Next First Last Back Forward 6
Theorem3.若一个链是一致的几儿何遍历转移概率以速度入，0<入<1收 敛的，以及相对于平稳分布是平方可积的，则有 vntn-Ef( →N0,1),n→o0 这个定理对构造平稳分布π的参数的置信区间提供了理论依据。 这些定理说明了： ·不可约十非周期+正常返一唯一的平稳分布. ·LLN:马氏链的实值函数的遍历均值以概率1收敛到极限分布下的均值. ·CLT:遍历均值作合适的变化依分布收敛到标准正态分布 因此，我们可以建立一个以π为平稳分布的马氏链，则在运行此链足够长时 间后，该马氏链就会达到平稳状态，此时马氏链的值就相当于从分布π中抽取 样本.从而可以根据遍历定理对积分进行估计.MCMC方法的优点在于： 1.适用于广泛且困难的问题： 2.问题维数的增加通常不会降低其收敛速度或者使其更复杂， Previous Next First Last Back Forward 7
Theorem 3. eòáÛ¥òóA¤H{(=£V«±Ñ›λ t , 0 < λ < 1¬ Ò), ±9 fÉÈu²©Ÿπ¥²êå», Kk √ n f¯n − Eπf(x) τ → N(0, 1), n → ∞ ˘á½nÈE²©ŸπÎÍò&´mJ¯ nÿù‚. ˘ ½n`² : • ÿå+ö±œ+~à=⇒ çò²©Ÿ. • LLN: ͺۢäºÍH{˛ä±V«1¬Ò4Å©Ÿe˛ä. • CLT: H{˛ää‹·Czù©Ÿ¬ÒIO©Ÿ. œd, ·Çå±Ô·òá±π販ŸÍºÛ, K3$1dÛv û m, Tͺۓ¨à ²G, dûͺÛä“Éul©Ÿπ•ƒ . l å±ä‚H{½nÈ»©?1O. MCMCê{`:3u: 1. ·^u2çÖ(JØK; 2. ØKëÍO\œ~ÿ¨¸$Ÿ¬ÒÑ›½ˆ¶ŸçE,. Previous Next First Last Back Forward 7
1.1.3 Markov Chain 离散的状态空间 假设{Xt,t≥0}为一齐次马氏链，转移概率为P,=P(Xt+1=Xt=), 平稳分布为π.则从相反的方向看此马氏链有 P(X:=jlX+1Xt+2...)=P(X=jlX+1), 即反方向来看，仍为马氏链(reversed markov chain,逆向马氏链).记P(t)为 反方向看马氏链从t十1时刻转移到t时刻的转移概率，则 P()= P(X:=X+1=)=PX+1=X:=)PX=) P(X+1=i P(】 开t+1() 这里的：表示马氏链在时刻的边际分布.从而逆向马氏链一般不再是齐次的， 如果当t→∞时其收敛到平稳分布，或者π0=π，则逆向马氏链就具有了平稳 Previous Next First Last Back Forward 8
1.1.3 Markov Chain l—Gòm b{Xt, t ≥ 0}èò‡gͺÛ, =£V«èPij = P(Xt+1 = j|Xt = i), ²©Ÿèπ. KlÉáêïwdͺÛk P(Xt = j|Xt+1, Xt+2, . . .) = P(Xt = j|Xt+1), =áêï5w, EèͺÛ(reversed markov chain, _ïͺÛ). PP ∗ ij (t)è áêïwͺÛlt + 1ûè=£tûè=£V«, K P ∗ ij (t) = P(Xt = j|Xt+1 = i) = P(Xt+1 = i|Xt = j)P(Xt = j) P(Xt+1 = i = Pjiπt(j) πt+1(i) ˘pπtL´ÍºÛ3tûè>S©Ÿ. l _ïͺÛòÑÿ2¥‡g. XJt → ∞ûŸ¬Ò²©Ÿ, ½ˆπ0 = π, K _ïͺۓ‰k ² Previous Next First Last Back Forward 8