R语言: Gibbs抽样实验

这篇文章并非原创，只是对下面PPT的总结和理解，作为入门性文字说明。
http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc.pdf

MH算法：Metropolis Hasting Algorithm
Gibbs抽样：Gibbs Sampling

上述两个算法是两个典型的马尔科夫链蒙特卡洛方法（Markov Chain Mento Carlo Method），作为随机抽样方法，普遍用于多元随机变量的抽样。通俗的说，通过迭代方法从目标随机概率分布中抽取样本。

关于蒙特卡洛方法和马尔科夫链，以及大数法则和中心极限定理，这边的证明略去。PPT中有详细解释。

1. Gibbs 抽样
Gibbs抽样的本质就是从已知的联合概率分布p(θ1, θ2,....|Y)（后验分布）推导出满条件分布，然后从满条件分布中抽取样本——p(θj|θ-j, Y)
这里有个问题：为什么可以从联合概率分布中求解满条件分布？
因为根据Hammersley-Clifford Theorem，任何联合概率可以由他的条件概率计算得来。
有了以上基本概念后，Gibbs抽样可以总结，先求取满条件分布，再对满条件分布抽烟。
满条件分布求取：
     a）先对后验分布进行整理，忽略相关的常数
     b）分别对参数θ1,θ2,.......进行条件分布求解
     c）正规化条件分布
Gibbs Sampler： 假设有三个参数θ1,θ2,θ3，后验分布为p(θ|Y)
     a）给参数 θ 初始值，记作 θ(0), θ取自于初始转移分布得来。
     b）对任意一个参数进行基于满条件分布的参数估计。为方便起见，我们从θ1开始，从p(θ1|θ2(0),θ3(0),Y)中得到θ1(1)
     c）同理，从p(θ2|θ1(1),θ3(0),Y)中得到θ2(1)
     d）从p(θ3|θ1(1),θ2(1),Y)中得到θ3(1)
     e）反复b,c,d步骤，直到求取到M个值之后停止。
     f）可以对样本进行burn-in和thinning操作
例子：
       一个核电厂有十个水泵，已知数据各水泵出故障次数，以及观察到出故障的时间
     

 
       y <- c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)
       t <- c(94, 16, 63, 126, 5, 31, 1, 1, 2, 10)
      

建立关于水泵故障次数的泊松分布，λi 为单位时间每个水泵的故障次数。似然函数有：
                   ∏i=1~10 Poisson(λi *ti)
假设Gamma(α, β) 是 λ 的先验分布。α=1.8. 所有的 λi都来自于这个分布。
假设Gamma(γ , δ) 是 β 的先验分布。γ=0.01 , δ=1
所以有：
p(λ, β|y, t) ∝  (∏i=1~10 Poisson(λi *ti)*Gamma(α, β))*Gamma(γ , δ)
经过整理有：
p(λ, β|y, t) ∝ ( ∏i=1~10 λi^(yi +α−1)*e^(−(ti +β)λi))*(β^(10α+γ−1)e^(−δβ))
通过上述联合分布，可以求得满条件分布:
p(λi |λ−i , β, y, t) ∝ λi^(yi +α−1)*e^(−(ti +β)λi)
p(β|λ, y, t)∝ β^(10α+γ−1)e^−β(δ+Sumi=1~10 λi)

gibbs<-function(n.sims,beta.start,alpha,gamma,delta,y,t,burnin=0,thin=1){
  beta.draws<-c()
  lambda.draws<-matrix(NA,nrow=n.sims,ncol=length(y))
  beta.cur<-beta.start
  lambda.update<-function(alpha,beta,y,t){
    rgamma(length(y),y+alpha,t+beta)
  }
  beta.update<-function(alpha,gamma,delta,lambda,y){
    rgamma(1,length(y)*alpha+gamma,delta+sum(lambda))
  }
  for(i in 1:n.sims){
    lambda.cur<-lambda.update(alpha=alpha,beta=beta.cur,y=y,t=t)
    beta.cur<- beta.update(alpha=alpha,gamma=gamma,delta=delta,lambda=lambda.cur,y=y)
    if(i>burnin&&(i-burnin)%%thin==0){
      lambda.draws[(i-burnin)%/%thin,]<-lambda.cur
      beta.draws[(i-burnin)/thin]<-beta.cur
    }
  }
  return(list(lambda.draws=lambda.draws,beta.draws=beta.draws))
}

其中lambda.update 方法和beta.update 方法是Gibbs Sampler过程。循环语句是对样本进行burn-in和thinning操作。

最后求取平均值，即所求参数λ, β的期望值。

posterior <- gibbs(n.sims = 10000, beta.start = 1, alpha = 1.8, gamma = 0.01, delta = 1, y = y, t = t)
colMeans(posterior$lambda.draws)
mean(posterior$beta.draws)