采样-Gibbs采样

MCMC蒙特卡洛马尔科夫链采样，非常重要的采样算法，而Gibbs算法也是MCMC种的一种，主要用于高维分布的采样。介绍MCMC的书籍有很多，https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/这是有关Gibbs采样matlab的一个实现，里面也介绍了gibbs采样的优缺点，缺点就是你得有高维分布它的条件分布的通项，如果它的条件分布不太好采样的话，可以使用拒绝采样来sample.Gibbs采样的算法我就不给了，因为网上有好多。

本文主要是在python上实现Gibbs采样，我们设定的是二维标准高斯分布的采样，其条件分布如下：

我画的图是动态图，代码如下：

 
 

   
   

    
    

   
   
   
   

    
    

     
     import matplotlib.pyplot 
     
     as plt
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     from scipy 
     
     import stats
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     import numpy 
     
     as np
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #为了方便操作，我们定义如下参数：
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #建立个字典吧
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #下面给的是标准正太分布，你可以修改参数
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     u1=
     
     0
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     u2=
     
     0
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     sigma_11=
     
     1
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     sigma_12=
     
     0
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     sigma_21=
     
     0
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     sigma_22=
     
     1
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     sample=[]
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     u=np.array([u1,u2],dtype=np.float32)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     sigma=np.array([[sigma_11,sigma_12],[sigma_21,sigma_22]],dtype=np.float32)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     def x1_given_x2(x2):
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         u_cond=u1+sigma_12/sigma_22*(x2-u2)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         sigma_cond=sigma_11-sigma_12/sigma_22*sigma_21
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     return stats.norm.rvs(loc=u_cond,scale=sigma_cond,size=
     
     1)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     def x2_given_x1(x1):
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         u_cond = u2 + sigma_12 / sigma_11 * (x1 - u1)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         sigma_cond = sigma_22 - sigma_12 / sigma_11 * sigma_12
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     return stats.norm.rvs(loc=u_cond, scale=sigma_cond, size=
     
     1)
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     def Gibbs_Sampling(initial,size):
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         sample.append(initial)
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     for i 
     
     in range(size):
    
    
   
   

   
   

    
    

   
   
   
   

    
    
        
     
     if i%
     
     2==
     
     0:
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
                 x1=x1_given_x2(sample[
     
     -1][
     
     1])
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
                 sample.append([x1,sample[
     
     -1][
     
     1]])
    
    
   
   

   
   

    
    

   
   
   
   

    
    
            
     
     continue
    
    
   
   

   
   

    
    

   
   
   
   

    
    
        
     
     else:
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
                 x2=x2_given_x1(sample[
     
     -1][
     
     0])
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
                 sample.append([sample[
     
     -1][
     
     0],x2])
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     return sample
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #初始样本点，我是自己定义的，你可以随机初始化
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     initial=np.array([
     
     5,
     
     8],dtype=np.float32)
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     b=Gibbs_Sampling(initial,
     
     151)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     b=np.array(b)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     for i 
     
     in range(len(b)):
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         plt.cla()
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         plt.plot(b[
     
     0:i+
     
     1,
     
     0],b[
     
     0:i+
     
     1,
     
     1])
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         plt.pause(
     
     0.03)
     
     #停顿时间，你可以调整
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     fig2=plt.figure()
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     ax=fig2.gca()
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     ax.scatter(b[:,
     
     0],b[:,
     
     1])
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     plt.show()
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   
 
 

采样过程的图片如下：

 

采样结束后得到的样本如下：