采样-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()





           

采样过程的图片如下：

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