Gibbs sampling & R

Gibbs sampling & R

分类： R 统计计算statistical computing 抽样、模拟、实验设计 2011-12-13 21:28  404人阅读  评论(0)  收藏  举报

维基百科的定义：http://en.wikipedia.org/wiki/Gibbs_sampling

我们假设一个随机的二元变量(x,y)，然后计算其中一个或全部的边缘分布p(x),p(y)。这种抽样思想是考虑条件分布p(xly),p(ylx)，比通过联合密度
p(xy)来计算求解简单,例如

首先，我们要进行初始化y0，然后通过条件概率P(xIy=y0)计算得出x0,然后再基于x0的条件分布p(ylx=x0)提取出一个更新后的y1。样本的更新过程是这样的:

,

重复该过程，可以生成一个GIBBS序列，其算法可抽象为：

例如：模拟一个(x,y)二元正态分布，该分布的一个边际分布为标准正态分布x~N(0,1)，x,y间的相关系数记为r。

我们利用Gibbs抽样进行模拟：

[html]  view plain copy

1. gibbs<-function (n, r)   
2. {  
3.         mat <- matrix(ncol = 2, nrow = n)  
4.         x <- 0  
5.         y <- 0  
6.         mat[1, ] <- c(x, y)  
7.         for (i in 2:n) {  
8.                 x <- rnorm(1, r * y, sqrt(1 - r^2))  
9.                 y <- rnorm(1, r * x, sqrt(1 - r^2))  
10.                 mat[i, ] <- c(x, y)  
11.         }  
12.         mat  
13. }  
14. bvn<-gibbs(10000,0.98)  

当然这个例子也可以直接模拟服从标准正态分布的边际分布x,再利用条件分布y|x进行模拟

[html]  view plain copy

1. rbvn<-function (n, r)   
2. {  
3.         x <- rnorm(n, 0, 1)  
4.         y <- rnorm(n, r * x, sqrt(1 - rho^2))  
5.         cbind(x, y)  
6. }  
7. bvn<-rbvn(10000,0.98)  

可以通过作图来看看两种模拟的效果

[html]  view plain copy

1. par(mfrow=c(3,2))  
2. plot(bvn,col=1:10000)  
3. plot(bvn,type="l")  
4. plot(ts(bvn[,1]))  
5. plot(ts(bvn[,2]))  
6. hist(bvn[,1],40)  
7. hist(bvn[,2],40)  
8. par(mfrow=c(1,1))  

结果如图

 

资料来源：

《基于网络的可适应Gibbs抽样方法研究》

http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/gibbs/gibbs.html

 

更多例子：

a.demo for one-way random effects model.

http://www.stat.sc.edu/~grego/courses/stat740/gibbs.1wayre.txt

b.simple Gibbs sampler class demonstration.

http://www.stat.sc.edu/~grego/courses/stat740/NWK1732.txt

 

更多资源：

《Introduction to Probability Simulation and Gibbs Sampling With R》

《Bayesian Computation With R》

《Introducing Monte Carlo Methods with R (Use R)》