中国科学技术大学：《实用统计软件》课程课件讲义（统计计算与软件）第九讲 Markov Chain Monte Carlo（二）马尔科夫蒙特卡罗方法

Lecture 9:Markov Chain Monte Carlo Methods (二) (马尔科夫蒙特卡罗方法) 张伟平 Monday 2nd November,2009

Lecture 9: Markov Chain Monte Carlo Methods () (ÍâÅÑAk¤ê{) ‹ï² Monday 2nd November, 2009

Contents 1 Markov Chain Monte Carlo Methods y 1.4 The Gibbs Sampler....................... 1.4.1 The Slice Gibbs Sampler................ 12 1.5 Monitoring Convergence······ 21 1.5.1.Convergence diagnostics plots·..·········· 21 1.5.2 Monte Carlo Error .... 22 1.5.3 The Gelman-Rubin Method ............. 25 l.6 VinBUGS Introduction·.···.. 33 l.6.1'Building Bayesian models in WinBUGS·.·.·.· 33 1.6.2 Model specification in WinBUGS..········· 36 1.6.3 Data and initial value specification。..·.·.··· 38 1.6.4 Compiling model and simulating values .. 46 Previous Next First Last Back Forward 1

Contents 1 Markov Chain Monte Carlo Methods 1 1.4 The Gibbs Sampler . . . . . . . . . . . . . . . . . . . . . . . 1 1.4.1 The Slice Gibbs Sampler . . . . . . . . . . . . . . . . 12 1.5 Monitoring Convergence . . . . . . . . . . . . . . . . . . . . 21 1.5.1 Convergence diagnostics plots . . . . . . . . . . . . . 21 1.5.2 Monte Carlo Error . . . . . . . . . . . . . . . . . . . 22 1.5.3 The Gelman-Rubin Method . . . . . . . . . . . . . . 25 1.6 WinBUGS Introduction . . . . . . . . . . . . . . . . . . . . 33 1.6.1 Building Bayesian models in WinBUGS . . . . . . . 33 1.6.2 Model specification in WinBUGS . . . . . . . . . . . 36 1.6.3 Data and initial value specification . . . . . . . . . . 38 1.6.4 Compiling model and simulating values . . . . . . . 46 Previous Next First Last Back Forward 1

Chapter 1 Markov Chain Monte Carlo Methods 1.4 The Gibbs Sampler Gibbs抽样是MH算法的一个特例，其经常用于目标分布是多元分布的场合. 假设所有的一元条件分布（每个分量对其他分量的条件分布）都是可以确定的， Gibbsi抽样使用这些一元条件分布进行抽样. 令X=(X1,·,Xd)为R中的随机变量，定义d-1维的随机变量 X-j=(X1,…,Xj-1,Xj+1,…,Xa, 并记XX-的条件密度为f(XX-).则Gibbs抽样是从这d个条件分布中 产生候选点.算法如下： 1.在t=0时，初始化X(0) Previous Next First Last Back Forward 1

Chapter 1 Markov Chain Monte Carlo Methods 1.4 The Gibbs Sampler Gibbs ƒ¥MHé{òáA~, Ÿ²~^u8I©Ÿ¥ı©Ÿ|‹. b§kò^á©Ÿ (zá©˛ÈŸ¶©˛^á©Ÿ)—¥å±(½, Gibbsƒ¶^˘ ò^á©Ÿ?1ƒ. -X = (X1, · · · , Xd)èRd•ëÅC˛, ½¬d − 1ëëÅC˛ X−j = (X1, · · · , Xj−1, Xj+1, · · · , Xd), øPXj |X−j^áó›èf(Xj |X−j ). KGibbsƒ¥l˘dá^á©Ÿ• )ˇ¿:. é{Xe: 1. 3t = 0û, –©zX(0); Previous Next First Last Back Forward 1

2.对t=1,2,·,T (a)令x1=X1(t-1). (b)对每个分量j=1,.,d, ()从f(Xz-)中产生候选点X() ()更新x=X() (c)今X(t)=(X(t),,X()(每个候选点都被接受) (d)增加t 注意在上述算法(b)步抽样中，各个分量依次被更新： x1(t)f(x1r2(t-1),…,xa(t-1): r2(t)~f(r2r1(t),x3(t-1),…,xd(t-1) xa(t)~f(zalzi(t),....xd_1(t)). Previous Next First Last Back Forward 2

2. Èt = 1, 2, · · · , T, (a) -x1 = X1(t − 1). (b) Èzᩲj = 1, . . . , d, (i) lf(Xj |x−j )•)ˇ¿:X∗ j (t). (ii) ç#xj = X∗ j (t). (c) -X(t) = (X∗ 1 (t), . . . , X∗ d (t))(záˇ¿:—…) (d) O\t 5ø3˛„é{(b)⁄ƒ•, àᩲùgç#: x1(t) ∼ f(x1|x2(t − 1), · · · , xd(t − 1)); x2(t) ∼ f(x2|x1(t), x3(t − 1), · · · , xd(t − 1)) . . . xd(t) ∼ f(xd|x1(t), · · · , xd−1(t)). Previous Next First Last Back Forward 2

从一元分布f(xz1(t),x2(t),…,xj-1(t),x5+1(t-1),…,x(t-1)中 抽样是比较容易的，因为f(xz-)xf(x),其中除了变量x,外，其他变量都是 常数 例1(Gibbs抽样：二元分布)使用Gibbs:抽样产生二元正态分布 N(41:2,ai,,P)的随机数 在二元正态场合，X1X2以及X2X1仍然服从正态分布，且易知 E[XlX2=x2l=h+p(2-2), 02 Var[X1|X2 =x2]=(1-p2)oi 类似可得X2X1的分布.因此 felx2)~N41+p(e2-42).(1-p)) 02 fe2lr)~N2+p2(e1-4,=1-p2) 因此，使用Gibbs算法如下 Previous Next First Last Back Forward 3

lò©Ÿf(xj |x1(t), x2(t), · · · , xj−1(t), xj+1(t − 1), · · · , xd(t − 1))• ƒ¥'N¥, œèf(xj |x−j ) ∝ f(x), Ÿ•ÿ C˛xj ,Ÿ¶C˛—¥ ~Í. ~ 1 (Gibbs ƒ: ©Ÿ) ¶^Gibbsƒ)©Ÿ N(µ1, µ2, σ2 1 , σ2 2 , ρ) ëÅÍ 3|‹, X1|X2±9X2|X1E,—l©Ÿ, Ö¥ E[X1|X2 = x2] = µ1 + ρ σ1 σ2 (x2 − µ2), V ar[X1|X2 = x2] = (1 − ρ 2 )σ 2 1 aqåX2|X1©Ÿ. œd f(x1|x2) ∼ N(µ1 + ρ σ1 σ2 (x2 − µ2), (1 − ρ 2 )σ 2 1 ) f(x2|x1) ∼ N(µ2 + ρ σ2 σ1 (x1 − µ1), = (1 − ρ 2 )σ 2 2 ) œd, ¶^Gibbsé{Xe Previous Next First Last Back Forward 3

1.令(x1,x2)=X(t-1)月 2.从f(x1E2)中产生候选点X(t): 3.更新x1=X(t): 4.从f(x21)中产生X5(t): 5.令X(t)=(X(),X(t): R代码如下 Cde #initialize constants and parameters N<-5000 #length of chain burn<-1000 #burn-in length X <-matrix(0,N,2) #the chain,a bivariate sample rho<--.75 #correlation mu1<-0 mu2<-2 sigma1 <-1 s1gma2<-.5 s1 <sqrt(1-rho2)*sigmal Previous Next First Last Back Forward 4

1. -(x1, x2) = X(t − 1); 2. lf(x1|x2)•)ˇ¿:X∗ 1 (t). 3. ç#x1 = X∗ 1 (t). 4. lf(x2|x1)•)X∗ 2 (t). 5. -X(t) = (X∗ 1 (t), X∗ 2 (t)). R ìËXe ↑Code #initialize constants and parameters N <- 5000 #length of chain burn<- 1000 #burn-in length X <- matrix(0, N, 2) #the chain, a bivariate sample rho <- -.75 #correlation mu1 <- 0 mu2 <- 2 sigma1 <- 1 sigma2 <- .5 s1 <- sqrt(1-rho^2)*sigma1 Previous Next First Last Back Forward 4

s2 <sqrt(1-rho 2)*sigma2 ####generate the chain##### X[1,]<-c(mu1,mu2) #initialize for (i in 2:N){ x2<-X[1-1,2] m1 <mu1 rho (x2 -mu2)*sigma1/sigma2 X[1,1]<-rnorm(1,m1,s1) x1<-X[1,1] m2 <mu2 rho (x1 mu1)*sigma2/sigma1 X[i,2]<rnorm(1,m2,s2) b <burn +1 x<-X[b:N,] Code 产生的链开始的1000个观测被丢弃掉，剩下的观测存在x中，对此样本计算均 值和协方差矩阵如下。各参数的样本估计离真值很近，散点图也显示出二元正 Previous Next First Last Back Forward 5

s2 <- sqrt(1-rho^2)*sigma2 ###### generate the chain ##### X[1, ] <- c(mu1, mu2) #initialize for (i in 2:N) { x2 <- X[i-1, 2] m1 <- mu1 + rho * (x2 - mu2) * sigma1/sigma2 X[i, 1] <- rnorm(1, m1, s1) x1 <- X[i, 1] m2 <- mu2 + rho * (x1 - mu1) * sigma2/sigma1 X[i, 2] <- rnorm(1, m2, s2) } b <- burn + 1 x <- X[b:N, ] ↓Code )Ûm©1000á*ˇøÔK, êe*ˇ3x•, Èd Oé˛ ä⁄ê› Xe. àÎÍOl˝äÈC, —:„èw´—  Previous Next First Last Back Forward 5

态所具有的球面对称性和负相关性特征 TCode compare sample statistics to parameters colMeans(x) cov(x) cor(x) plot(x,main="",cex=.5,xlab=bquote(X[1]), ylab=bquote(X[2]),ylim=range(x[,2])) Code 例2（贝叶斯分析例子：身体温度数据）考虑Mackowiak et al. (1992)的数据，该数据记录了130个人的身体温度（华氏），性别和每分钟的 心跳.实验的目的是检验Carl Wunderlich的观点-健康成年人的体温平 均为37C(=98.6°F). Previous Next First Last Back Forward 6

§‰k•°È°5⁄KÉ'5A. ↑Code # compare sample statistics to parameters colMeans(x) cov(x) cor(x) plot(x, main="", cex=.5, xlab=bquote(X[1]), ylab=bquote(X[2]), ylim=range(x[,2])) ↓Code ~ 2 (ìd©¤~f: Nߛ͂) ƒ Mackowiak et al. (1992)Í‚, TÍ‚P¹ 130á<Nß›(uº), 5O⁄z©® %a. ¢8¥uCarl Wunderlich*:– Ëx§c<Nß² ˛è37oC(=98.6oF). Previous Next First Last Back Forward 6

记温度为，i=1,·,n,并假设正态模型 V(4,a2) 以及取先验分布为 N(uo,ao),2~IG(ao,bo) 则此时我们的目标分布为4，σ的后验分布 f(4,o21g)xf(4,a2)π(4)r(a2) 为使用Gibbs抽样算法，我们必须计算f(叫a2,)与f(a2|4,y).经过计算我们 得到 哈 4o,g≈Ng+=o,wnw=2n+品 n ou.y1Ga+号o+∑-P. n 2 i=1 Previous Next First Last Back Forward 7

Pß›èyi, i = 1, · · · , n, øb. yi ∼ N(µ, σ2 ) ±9k©Ÿè µ ∼ N(µ0, σ0), σ2 ∼ IG(a0, b0) Kdû·Ç8I©Ÿèµ, σ2￾©Ÿ f(µ, σ2 |y) ∝ f(y|µ, σ2 )π(µ)π(σ 2 ) è¶^Gibbsƒé{, ·Ç7LOéf(µ|σ 2 , y)Üf(σ 2 |µ, y). ²LOé·Ç  µ|σ, y ∼ N(ωy¯ + (1 − ω)µ0, ω σ 2 n ), ω = σ 2 0 σ2/n + σ 2 0 . σ 2 |µ, y ∼ IG(a0 + n 2 , b0 + 1 2 Xn i=1 (yi − µ) 2 ). Previous Next First Last Back Forward 7

使用此结果，Gibbs抽样算法如下： 对t=1,…,T, 1.令μ=ht-1),g=at-1) 2.计第w=7m=u呵+(1-ω)o和s2=w民 3.从N(m,s2)中产生μ. 4.令u(8=μ 5.计算a=a0+受，b=b0+∑=1(h-4)2 6.从G(a,b)中产生T 7.令g2=1/r以及(0=. 从而R代码如下： 下coda bodytemp<-read.table("bodytemp.txt",header=T) y<-bodytempStemp bary<-mean(y);n<-length(y) Iterations<-3500 mu0<-0;s0<-100;a0<-0.001;b0<-0.001 Previous Next First Last Back Forward 8

¶^d(J, Gibbsƒé{Xe: Èt = 1, · · · , T, 1. -µ = µ (t−1), σ = σ (t−1) . 2. Oéω = σ 2 0 σ2/n+σ2 0 , m = ωy¯ + (1 − ω)µ0 ⁄s 2 = ω σ 2 n . 3. lN(m, s2 )•)µ. 4. -µ (t) = µ. 5. Oéa = a0 + n 2 , b = b0 + 1 2 Pn i=1(yi − µ) 2 . 6. lG(a, b)•)τ. 7. -σ 2 = 1/τ±9σ (t) = σ. l RìËXe: ↑Code bodytemp<-read.table("bodytemp.txt",header=T) y<-bodytemp$temp bary<-mean(y); n<-length(y) Iterations<-3500 mu0<-0; s0<-100; a0<-0.001; b0<-0.001 Previous Next First Last Back Forward 8