MCMC采样_MCMC认证

大家好，又见面了，我是你们的朋友全栈君。

MCMC(一)蒙特卡罗方法

MCMC(二)马尔科夫链

MCMC(三)MCMC采样和M-H采样

　　　　MCMC(四)Gibbs采样

　　　　在MCMC(三)MCMC采样和M-H采样中，我们讲到了M-H采样已经可以很好的解决蒙特卡罗方法需要的任意概率分布的样本集的问题。但是M-H采样有两个缺点：一是需要计算接受率，在高维时计算量大。并且由于接受率的原因导致算法收敛时间变长。二是有些高维数据，特征的条件概率分布好求，但是特征的联合分布不好求。因此需要一个好的方法来改进M-H采样，这就是我们下面讲到的Gibbs采样。

1. 重新寻找合适的细致平稳条件

满足：

　　　　则称概率分布\pi(x)是状态转移矩阵P的平稳分布。

　　　　在M-H采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。

，容易发现下面两式成立：

　　　　由于两式的右边相等，因此我们有：

　　　　也就是：

,我们可以得到：

：

，满足细致平稳条件时：

　　　　于是这个二维空间上的马氏链将收敛到平稳分布 \pi(x,y)

2. 二维Gibbs采样

　　　　利用上一节找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样需要两个维度之间的条件概率。具体过程如下：

　　　　1）输入平稳分布\pi(x_1,x_2)，设定状态转移次数阈值n_1，需要的样本个数n_2

　　　　2）随机初始化初始状态值x_1^{(0)}和x_2^{(0)}

　　　　3）for t = 0 to n_1 +n_2-1:

　　　　　　a) 从条件概率分布P(x_2|x_1^{(t)})中采样得到样本x_2^{t+1}

　　　　　　b) 从条件概率分布P(x_1|x_2^{(t+1)})中采样得到样本x_1^{t+1}

　　　　样本集\{(x_1^{(n_1)}, x_2^{(n_1)}), (x_1^{(n_1+1)}, x_2^{(n_1+1)}), …, (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)})\}即为我们需要的平稳分布对应的样本集。

　　　　整个采样过程中，我们通过轮换坐标轴，采样的过程为：

　　　　用下图可以很直观的看出，采样是在两个坐标轴上不停的轮换的。当然，坐标轴轮换不是必须的，我们也可以每次随机选择一个坐标轴进行采样。不过常用的Gibbs采样的实现都是基于坐标轴轮换的。

MCMC采样_MCMC认证

3. 多维Gibbs采样

　　　　上面的这个算法推广到多维的时候也是成立的。比如一个n维的概率分布\pi(x_1,x_2,…x_n)，我们可以通过在n个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴x_i上的转移，马尔科夫链的状态转移概率为P(x_i|x_1,x_2,…,x_{i-1},x_{i+1},…,x_n)，即固定n-1个坐标轴，在某一个坐标轴上移动。

　　　　具体的算法过程如下：

　　　　1）输入平稳分布\pi(x_1,x_2，…,x_n)或者对应的所有特征的条件概率分布，设定状态转移次数阈值n_1，需要的样本个数n_2

　　　　2）随机初始化初始状态值(x_1^{(0)},x_2^{(0)},…,x_n^{(0)}

　　　　3）for t = 0 to n_1 +n_2-1:

　　　　　　a) 从条件概率分布P(x_1|x_2^{(t)}, x_3^{(t)},…,x_n^{(t)})中采样得到样本x_1^{t+1}

　　　　　　b) 从条件概率分布P(x_2|x_1^{(t+1)}, x_3^{(t)}, x_4^{(t)},…,x_n^{(t)})中采样得到样本x_2^{t+1}

　　　　　　c)…

　　　　　　d) 从条件概率分布P(x_j|x_1^{(t+1)}, x_2^{(t+1)},…, x_{j-1}^{(t+1)},x_{j+1}^{(t)}…,x_n^{(t)})中采样得到样本x_j^{t+1}

　　　　　　e)…

　　　　　　f) 从条件概率分布P(x_n|x_1^{(t+1)}, x_2^{(t+1)},…,x_{n-1}^{(t+1)})中采样得到样本x_n^{t+1}

　　　　样本集\{(x_1^{(n_1)}, x_2^{(n_1)},…, x_n^{(n_1)}), …, (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)},…,x_n^{(n_1+n_2-1)})\}即为我们需要的平稳分布对应的样本集。

　　　　整个采样过程和Lasso回归的坐标轴下降法算法非常类似，只不过Lasso回归是固定n-1个特征，对某一个特征求极值。而Gibbs采样是固定n-1个特征在某一个特征采样。

　　　　同样的，轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。

4. 二维Gibbs采样实例

　　　　这里给出一个Gibbs采样的例子。完整代码参见我的github: https://github.com/ljpzzz/machinelearning/blob/master/mathematics/mcmc_3_4.ipynb

,其中：

　　　　而采样过程中的需要的状态转移条件分布为：

　　　　具体的代码如下：

from mpl_toolkits.mplot3d import Axes3D from scipy.stats import multivariate_normal samplesource = multivariate_normal(mean=[5,-1], cov=[[1,1],[1,4]]) def p_ygivenx(x, m1, m2, s1, s2): return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt((1 - rho ** 2) * (s2**2)))) def p_xgiveny(y, m1, m2, s1, s2): return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt((1 - rho ** 2) * (s1**2)))) N = 5000 K = 20 x_res = [] y_res = [] z_res = [] m1 = 5 m2 = -1 s1 = 1 s2 = 2 rho = 0.5 y = m2 for i in xrange(N): for j in xrange(K): x = p_xgiveny(y, m1, m2, s1, s2) y = p_ygivenx(x, m1, m2, s1, s2) z = samplesource.pdf([x,y]) x_res.append(x) y_res.append(y) z_res.append(z) num_bins = 50 plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5) plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5) plt.title('Histogram') plt.show()

　　　　输出的两个特征各自的分布如下：

MCMC采样_MCMC认证

　　　　然后我们看看样本集生成的二维正态分布，代码如下：

fig = plt.figure() ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20) ax.scatter(x_res, y_res, z_res,marker='o') plt.show()

　　　　输出的正态分布图如下：

MCMC采样_MCMC认证

5. Gibbs采样小结

　　　　由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。

　　　　有了Gibbs采样来获取概率分布的样本集，有了蒙特卡罗方法来用样本集模拟求和，他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。MCMC系列就在这里结束吧。

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/167285.html原文链接：https://javaforall.cn