很多时候,我们希望从后验分布中采样,用以估算后验分布的参数或其他用途。本文主要介绍后验分布的gipps采样方法。
假定我们有K个变量(U1,U2,...,UK)。希望从这K个变量的联合分布中采样。假定这是一个很困难的问题,但是从条件分布P(Ui|U(-i))中采样很方便。Gibbs采样就是不断的从条件分布中对各个Ui(i∈(1,2,...,K))分量采样,然后用采样的值去替换原有的变量Ui,直到该过程稳定下来。这里的i 取值范围是1-K. 条件分布P(Ui|U(-i))中的U(-i)是i取值1-K,但是不包括i。比如K=3,当i=2时,条件分布写成P(U2|U1,U3)。
具体举个例子,假定我们有三个变量U1,U2,U3,其联合分布为Pr(U1,U2,U3)。在步骤t,我们有值U1(t),U2(t),U3(t)。在t+1步,我们首先用新值U1(t+1)替换U1(t),U1(t+1)是基于条件分布P(U1|U2(t),U3(t))采样得到的。然后,我们根据P(U2|U1(t+1),U3(t))采样替换U2(t)。(注意,这里的条件分布已用U1(t+1)替换掉了U1(t))。接下来,我们用P(U3|U1(t+1),U2(t+1))采样替换U3(t).依次循环,直到收敛。采用的算法流程可描述如下:
Gipps采样的收敛性本文不解释。需要注意的是,我们并不需要知道P(Ui|U(-i))的确切形式,只要能从中采样就可以了。
书上给了个Gibbs采样和EM算法的关系的例子。还是以EM那一节混合高斯的例子来介绍。如下:
Gibbs采样中把隐变量看成是一个参数,于是混合高斯的参数为(theta,Zm)。为了简单起见,我们假定theta中的方差和混合比例pi都是已知的,未知的参数是混合高斯的均值(mu1和mu2)和隐变量Zm。Gipps采样估参的过程如算法8.4:
从中可以看到步骤2的a,b两步和EM算法非常类似。下面我们贴出书上图8.8的代码实现。
mu1 <-4;mu2 <- 1;pi <- 0.55;sigma1 <- 0.9;sigma2 <- 0.9;
y=c(-0.39, 0.12,0.94 ,1.67, 1.76, 2.44 ,3.72, 4.28, 4.92 ,5.53,0.06 ,0.48, 1.01, 1.68, 1.80 ,3.25, 4.12, 4.60, 5.28 ,6.22 )
m1 <- c()
m2 <- c()
for (i in 1:1000)
{
fhtheta1 <- dnorm(y,mean=mu1,sd=sigma1)
fhtheta2 <- dnorm(y,mean=mu2,sd=sigma2)
delta <- pi*fhtheta2/((1-pi)*fhtheta1+pi*fhtheta2)
delta <- ifelse(delta>0.5,1,0)
mu1Hat<- sum((1-delta)*y)/sum(1-delta)
mu2Hat <- sum(delta*y)/sum(delta)
mu1 <- rnorm(1,mu1Hat,sigma1)
mu2 <- rnorm(1,mu2Hat,sigma2)
m1 <- c(m1,mu1)
m2 <- c(m2,mu2)
}
plot(1:1000,m1,type="l",col="orange",ylim=c(-1,8))
lines(m2,col="blue")
但是这段代码执行之后的结果很奇怪, 有的时候能和书上的一致,有的时候不行,经常是mu1和mu2的估值互换,如下图所示,我检查了好几遍代码,始终没有找到问题所在。希望以后有人能够告诉我。