1. 重新寻找合适的细致平稳条件
在上一篇中,我们讲到了细致平稳条件:如果非周期马尔科夫链的状态转移矩阵P和概率分布π(x)对于所有的i,j满足:
则称概率分布π(x)是状态转移矩阵P的平稳分布。
在M-H采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。
从二维的数据分布开始,假设π(x1,x2)是一个二维联合数据分布,观察第一个特征维度相同的两个点A(x(1)1,x(1)2)和B(x(1)1,x(2)2),容易发现下面两式成立:
由于两式的右边相等,因此我们有:
也就是:
观察上式再观察细致平稳条件的公式,我们发现在x1=x(1)1这条直线上,如果用条件概率分布π(x2|x(1)1)作为马尔科夫链的状态转移概率,则任意两个点之间的转移满足细致平稳条件!这真是一个开心的发现,同样的道理,在在x2=x(1)2这条直线上,如果用条件概率分布π(x1|x(1)2)作为马尔科夫链的状态转移概率,则任意两个点之间的转移也满足细致平稳条件。那是因为假如有一点C(x(2)1,x(1)2),我们可以得到:
基于上面的发现,我们可以这样构造分布π(x1,x2)的马尔可夫链对应的状态转移矩阵P:
有了上面这个状态转移矩阵,我们很容易验证平面上的任意两点E,F,满足细致平稳条件:
2. 二维Gibbs采样
利用上一节找到的状态转移矩阵,我们就得到了二维Gibbs采样,这个采样需要两个维度之间的条件概率。具体过程如下:
1)输入平稳分布π(x1,x2),设定状态转移次数阈值n1,需要的样本个数n2
2)随机初始化初始状态值x(1)1和x(1)2
3)for t=0 to n1+n2−1:
a) 从条件概率分布P(x2|x(t)1)中采样得到样本xt+12
b) 从条件概率分布P(x1|x(t+1)2)中采样得到样本xt+11
样本集{(x(n1)1,x(n1)2),(x(n1+1)1,x(n1+1)2),...,(x(n1+n2−1)1,x(n1+n2−1)2)}即为我们需要的平稳分布对应的样本集。
整个采样过程中,我们通过轮换坐标轴,采样的过程为:
用下图可以很直观的看出,采样是在两个坐标轴上不停的轮换的。当然,坐标轴轮换不是必须的,我们也可以每次随机选择一个坐标轴进行采样。不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
3. 多维Gibbs采样
上面的这个算法推广到多维的时候也是成立的。比如一个n维的概率分布π(x1,x2,...xn),我们可以通过在n个坐标轴上轮换采样,来得到新的样本。对于轮换到的任意一个坐标轴xi上的转移,马尔科夫链的状态转移概率为P(xi|x1,x2,...,xi−1,xi+1,...,xn),即固定n−1个坐标轴,在某一个坐标轴上移动。
具体的算法过程如下:
1)输入平稳分布π(x1,x2,...,xn)或者对应的所有特征的条件概率分布,设定状态转移次数阈值n1,需要的样本个数n2
2)随机初始化初始状态值(x(1)1,x(1)2,...,x(1)n
3)for t=0 to n1+n2−1:
a) 从条件概率分布P(x1|x(t)2,x(t)3,...,x(t)n)中采样得到样本xt+11
b) 从条件概率分布P(x2|x(t+1)1,x(t)3,x(t)4,...,x(t)n)中采样得到样本xt+12
c)...
d) 从条件概率分布P(xj|x(t+1)1,x(t+1)2,...,x(t+1)j−1,x(t)j+1...,x(t)n)中采样得到样本xt+1j
e)...
f) 从条件概率分布P(xn|x(t+1)1,x(t+1)2,...,x(t+1)n−1)中采样得到样本xt+1n
样本集{(x(n1)1,x(n1)2,...,x(n1)n),...,(x(n1+n2−1)1,x(n1+n2−1)2,...,x(n1+n2−1)n)}即为我们需要的平稳分布对应的样本集。
整个采样过程和Lasso回归的坐标轴下降法算法非常类似,只不过Lasso回归是固定n−1个特征,对某一个特征求极值。而Gibbs采样是固定n−1个特征在某一个特征采样。
同样的,轮换坐标轴不是必须的,我们可以随机选择某一个坐标轴进行状态转移,只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
4. 二维Gibbs采样实例
这里给出一个Gibbs采样的例子。假设我们要采样的是一个二维正态分布Norm(μ,Σ),其中:
而采样过程中的需要的状态转移条件分布为:
具体的代码如下:
from mpl_toolkits.mplot3d import Axes3D from scipy.stats import multivariate_normal samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]]) def p_ygivenx(x, m1, m2, s1, s2): return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2)) def p_xgiveny(y, m1, m2, s1, s2): return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1)) 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()
输出的两个特征各自的分布如下:
然后我们看看样本集生成的二维正态分布,代码如下:
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()
输出的正态分布图如下:
5. Gibbs采样小结
由于Gibbs采样在高维特征时的优势,目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的,同时Gibbs采样要求数据至少有两个维度,一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。
有了Gibbs采样来获取概率分布的样本集,有了蒙特卡罗方法来用样本集模拟求和,他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。MCMC系列就在这里结束吧。
本文转自刘建平Pinard博客园博客,原文链接:http://www.cnblogs.com/pinard/p/6645766.html,如需转载请自行联系原作者