MCMC采样和M-H采样
第三章 矩阵和线性代数

机器学习原理

Gibbs采样

M-H采样有两个缺点：一是需要计算接受率，在高维时计算量大。并且由于接受率的原因导致算法收敛时间变长。二是有些高维数据，特征的条件概率分布好求，但是特征的联合分布不好求。因此需要一个好的方法来改进M-H采样，这就是我们下面讲到的Gibbs采样。

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

如果非周期马尔科夫链的状态转移矩阵P和概率分布$$\pi(x)$$对于所有的i,j满足：$$\pi(i)P(i,j) = \pi(j)P(j,i)$$

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

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

从二维的数据分布开始，假设$$\pi(x_1,x_2)$$是一个二维联合数据分布，观察第一个特征维度相同的两个点$$A(x_1{(1)},x_2{(1)})$$和$$B(x_1{(1)},x_2{(2)})$$，容易发现下面两式成立：$$\pi(x_1{(1)},x_2{(1)}) \pi(x_2{(2)} | x_1{(1)}) = \pi(x_1{(1)})\pi(x_2{(1)}|x_1{(1)}) \pi(x_2{(2)} | x_1{(1)})\pi(x_1{(1)},x_2{(2)}) \pi(x_2{(1)} | x_1{(1)}) = \pi(x_1{(1)}) \pi(x_2{(2)} | x_1{(1)})\pi(x_2{(1)}|x_1{(1)})$$

由于两式的右边相等，因此我们有：$$\pi(x_1{(1)},x_2{(1)}) \pi(x_2{(2)} | x_1{(1)}) = \pi(x_1{(1)},x_2{(2)}) \pi(x_2{(1)} | x_1{(1)})$$

也就是：$$\pi(A) \pi(x_2{(2)} | x_1{(1)}) = \pi(B) \pi(x_2{(1)} | x_1{(1)})$$

观察上式再观察细致平稳条件的公式，我们发现在$$x_1 = x_1{(1)}$$这条直线上，如果用条件概率分布$$\pi(x_2| x_1{(1)})$$作为马尔科夫链的状态转移概率，则任意两个点之间的转移满足细致平稳条件！这真是一个开心的发现，同样的道理，在在$$x_2 = x_2{(1)}$$这条直线上，如果用条件概率分布$$\pi(x_1| x_2{(1)})$$作为马尔科夫链的状态转移概率，则任意两个点之间的转移也满足细致平稳条件。那是因为假如有一点$$C(x_1{(2)},x_2{(1)})$$,我们可以得到：$$\pi(A) \pi(x_1{(2)} | x_2{(1)}) = \pi(C) \pi(x_1{(1)} | x_2{(1)})$$

基于上面的发现，我们可以这样构造分布$$\pi(x_1,x_2)$$的马尔可夫链对应的状态转移矩阵

$$P(A \to B) = \pi(x_2{(B)}|x_1{(1)});; if; x_1{(A)} = x_1{(B)} =x_1^{(1)}$$

$$P(A \to C) = \pi(x_1{(C)}|x_2{(1)});; if; x_2{(A)} = x_2{(C)} =x_2^{(1)}$$

$$P(A \to D) = 0;; else$$

有了上面这个状态转移矩阵，我们很容易验证平面上的任意两点E,F，满足细致平稳条件：$$\pi(E)P(E \to F) = \pi(F)P(F \to E)$$

2. 二维Gibbs采样

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

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

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

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)})}$$即为我们需要的平稳分布对应的样本集。

整个采样过程中，我们通过轮换坐标轴，采样的过程为：$$(x_1{(1)}, x_2{(1)}) \to (x_1{(1)}, x_2{(2)}) \to (x_1{(2)}, x_2{(2)}) \to ... \to (x_1{(n_1+n_2-1)}, x_2{(n_1+n_2-1)})$$

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

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{(1)},x_2{(1)},...,x_n^{(1)}$$

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采样的例子。假设我们要采样的是一个二维正态分布$$Norm(\mu,\Sigma)$$,其中：

$$\mu = (\mu_1,\mu_2) = (5,-1)$$

$$\Sigma = \left( \begin{array}{ccc} \sigma_12&\rho\sigma_1\sigma_2 \ \rho\sigma_1\sigma_2 &\sigma_22 \end{array} \right) = \left( \begin{array}{ccc} 1&1 \ 1&4 \end{array} \right)$$

而采样过程中的需要的状态转移条件分布为：$$P(x_1|x_2) = Norm\left ( \mu _1+\rho \sigma_1/\sigma_2 \left ( x _2-\mu _2 \right ), (1-\rho 2)\sigma_12 \right )P(x_2|x_1) = Norm\left ( \mu _2+\rho \sigma_2/\sigma_1 \left ( x _1-\mu _1 \right ), (1-\rho 2)\sigma_22 \right )$$

具体的代码如下：

import random
import math
import matplotlib.pyplot as plt
%matplotlib inline
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 * math.sqrt(s2 / s1) * (x - m1), math.sqrt(1 - rho ** 2) * s2))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * math.sqrt(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 = 4

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('Gibbs Sampler')
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算法在大数据时代高维数据模拟求和时的作用。

书籍推荐