ESL-chapter8-gibbs采样

     很多时候，我们希望从后验分布中采样，用以估算后验分布的参数或其他用途。本文主要介绍后验分布的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的估值互换，如下图所示，我检查了好几遍代码，始终没有找到问题所在。希望以后有人能够告诉我。