theta <-matrix(nrow=Iterations,ncol=2) cur.mu<-0;cur.tau<-2;cur.s<-sqrt(1/cur.tau) for (t in 1:Iterations){ <-s0^2/(cur.s2/n+s0^2) m <w*bary +(1-v)*muo s <sqrt(w/n )cur.s cur.mu <-rnorm(1,m,s a<-a0+0.5n b <b0 0.5 sum((y-cur.mu)2 cur.tau <-rgamma(1,a,b cur.s <sqrt(1/cur.tau) theta[t,]<-c(cur.mu,cur.s) mcmc.output<-theta apply(mcmc.output [-(1:1000),],2,mean) #compare to true value:98.25,0.542 apply(mcmc.output [-(1:1000),],2,sd) #compare to true value:0.06456,0.06826 Code Previous Next First Last Back Forward 9
theta <- matrix(nrow=Iterations, ncol=2) cur.mu<-0; cur.tau<-2; cur.s<-sqrt(1/cur.tau) for (t in 1:Iterations){ w<- s0^2/( cur.s^2/n+ s0^2 ) m <- w*bary + (1-w)*mu0 s <- sqrt( w/n ) * cur.s cur.mu <- rnorm( 1, m, s ) a <- a0 + 0.5*n b <- b0 + 0.5 * sum( (y-cur.mu)^2 ) cur.tau <- rgamma( 1, a, b ) cur.s <- sqrt(1/cur.tau) theta[t,]<-c( cur.mu, cur.s) } mcmc.output<-theta apply(mcmc.output[-(1:1000),],2,mean) #compare to true value: 98.25, 0.542 apply(mcmc.output[-(1:1000),],2,sd) #compare to true value: 0.06456, 0.06826 ↓Code Previous Next First Last Back Forward 9
对链的诊断图 par(mfrow=c(3,2),xaxs='r',yaxs='r',bty='1',cex=0.8) 1ter<-1500 burnin<-500 index<-1:iter index2<-(burnin+1):iter plot(index,theta[index,1],type='1',ylab='Values of mu', xlab='Iterations',main='(a)Trace Plot of mu') plot(index,theta[index,2],type='1',ylab='Values of sigma', xlab='Iterations',main='(b)Trace Plot of sigma') ergtheta0<-erg.mean(theta[index,1] ergtheta02<-erg.mean(theta[index2,1] ylims0<-range(c(ergtheta0,ergtheta02)) ergthetal<-erg.mean(theta[index,2] ergtheta12<-erg.mean(theta[index2,2] ylims1<-range(c(ergthetal,ergtheta12)) Previous Next First Last Back Forward 10
ÈÛ‰„ ↑Code par( mfrow=c(3,2), xaxs=’r’, yaxs=’r’, bty=’l’ , cex=0.8) iter<-1500 burnin<-500 index<-1:iter index2<-(burnin+1):iter plot(index, theta[index,1], type=’l’, ylab=’Values of mu’, xlab=’Iterations’, main=’(a) Trace Plot of mu’) plot(index, theta[index,2], type=’l’, ylab=’Values of sigma’, xlab=’Iterations’, main=’(b) Trace Plot of sigma’) ergtheta0<-erg.mean( theta[index,1] ) ergtheta02<-erg.mean( theta[index2,1] ) ylims0<-range( c(ergtheta0,ergtheta02) ) ergtheta1<-erg.mean( theta[index,2] ) ergtheta12<-erg.mean( theta[index2,2] ) ylims1<-range( c(ergtheta1,ergtheta12) ) Previous Next First Last Back Forward 10
step<-10 index3<-seq(1,iter,step) index4<-seq(burnin+1,iter,step) plot(index3 ergthetao[index3],type='1',ylab='Values of mu', xlab='Iterations',main='(c)Ergodic Mean Plot of mu',ylim=ylims0) lines(index4,ergtheta02[index4-burnin],col=2,lty=2) plot(index3,ergthetal[index3],type='1',ylab='Values of sigma', xlab='Iterations',main='(d)Ergodic Mean Plot of sigma',ylim=ylims1) lines(index4,ergtheta12[index4-burnin],col=2,lty=2) acf(theta[index2,1],main='Autocorrelations Plot for mu') acf(theta[index2,2],main='Autocorrelations Plot for sigma') Code Previous Next First Last Back Forward 11
step<-10 index3<-seq(1,iter,step) index4<-seq(burnin+1,iter,step) plot(index3 , ergtheta0[index3], type=’l’, ylab=’Values of mu’, xlab=’Iterations’, main=’(c) Ergodic Mean Plot of mu’, ylim=ylims0) lines(index4, ergtheta02[index4-burnin], col=2, lty=2) plot(index3, ergtheta1[index3], type=’l’, ylab=’Values of sigma’, xlab=’Iterations’, main=’(d) Ergodic Mean Plot of sigma’, ylim=ylims1) lines(index4, ergtheta12[index4-burnin], col=2, lty=2) acf(theta[index2,1], main=’Autocorrelations Plot for mu’) acf(theta[index2,2], main=’Autocorrelations Plot for sigma’) ↓Code Previous Next First Last Back Forward 11
1.4.1 The Slice Gibbs Sampler 切片Gibbs抽样(Slice Gibbs Sampler)本质上是基于Gibbs抽样的.其主要用 于当完全的条件分布没有简单或者方便的形式情形.这个方法通过添加一些 辅助变量把参数空间扩大，但保持感兴趣的边际分布不变，而把所有的条件分 布转变成标准形式.然后可以使用标准的Gbbs抽样方法. 切片Gibbs抽样的想法如下.考虑目标分布g(x),其很难进行抽样.我们引 入一个新的辅助变量u,其条件分布为f(uz).则联合分布为 f(u,x)=f(ulr)g(x) 而x的边际分布等于目标分布g(x).因此我们可以使用Gibbs算法从联合分布 f(u,x),以及边际分布f(u)和f(x)=g(x)中产生随机数： 1.产生uf(r) 2.产生x~f(r)g(r). 由于f(ur)在上述计算中出现两次，因此其的选取要使得从分布f(ux)和f(ux)g(x) Previous Next First Last Back Forward 12
1.4.1 The Slice Gibbs Sampler É°Gibbsƒ(Slice Gibbs Sampler) ü˛¥ƒuGibbsƒ. ŸÃá^ u ^á©Ÿvk{¸½ˆêB/™ú/. ˘áê{œLV\ò 9œC˛ rÎÍòm*å, ±a,>S©ŸÿC, r§k^á© Ÿ=C§IO/™. ,屶^IOGibbsƒê{. É°Gibbsƒé{Xe. ƒ8I©Ÿg(x), ŸÈJ?1ƒ. ·Ç⁄ \òá#9œC˛u, Ÿ^á©Ÿèf(u|x). KÈ‹©Ÿè f(u, x) = f(u|x)g(x) x>S©Ÿ u8I©Ÿg(x). œd·Ç屶^Gibbsé{lÈ‹©Ÿ f(u, x), ±9>S©Ÿf(u)⁄f(x) = g(x)•)ëÅÍ: 1. )u ∼ f(u|x). 2. )x ∼ f(u|x)g(x). duf(u|x)3˛„Oé•—y¸g, œdŸ¿á¶l©Ÿf(u|x)⁄f(u|x)g(x) Previous Next First Last Back Forward 12
中很方便的抽样，常用的一个选择是均匀分布U(O,g(x),此时 f(u,x)=- a96@10<u<96)=10<u<g》, f)= I(0<u<g(x))du=g(x). 因此，此时Gbbs抽样算法如下 1.产生因~U(0,9(x(-1)月 2.产生x~U(x:0≤u(≤9(x) 对贝叶斯分析来说，经常选择u~U(0,f(),联合分布为 fo,u)）xΠI0≤u:≤flof(o) 从而Gibbs算法如下 1.令0=8(t-10. 2.对i=1,,n,产生u9心U(0,f(l0) 3.对j=1,·,d,更新0，~f(0)Π11(0≤u9≤fg:l) Previous Next First Last Back Forward 13
•ÈêBƒ, ~^òá¿J¥˛!©ŸU(0, g(x)), dû f(u, x) = 1 g(x) g(x)I(0 < u < g(x)) = I(0 < u < g(x)), f(x) = Z I(0 < u < g(x))du = g(x). œd, dûGibbsƒé{Xe 1. )u (t) ∼ U(0, g(x (t−1))); 2. )x (t) ∼ U(x : 0 ≤ u (t) ≤ g(x)). Èìd©¤5`, ²~¿Ju ∼ U(0, f(y|θ)), È‹©Ÿè f(θ, u|y) ∝ {Yn i=1 I(0 ≤ ui ≤ f(yi|θ))}f(θ) l Gibbsé{Xe 1. -θ = θ (t−1) . 2. Èi = 1, . . . , n, )u (t) i ∼ U(0, f(yi|θ)). 3. Èj = 1, · · · , d, ç#θj ∼ f(θj ) Qn i=1 I(0 ≤ u (t) i ≤ f(yi|θ)). Previous Next First Last Back Forward 13