在《MCMC(二)：MCMC采样和M-H采样》中，我们讲到了M-H采样已经可以很好的解决蒙特卡罗方法需要的任意概率分布的样本集的问题。
但是M-H采样有两个缺点：一是需要计算接受率，在高维时计算量大。并且由于接受率的原因导致算法收敛时间变长。二是有些高维数据，特征的条件概率分布好求，但是特征的联合分布不好求。因此需要一个好的方法来改进M-H采样，这就是我们下面讲到的Gibbs采样。
本篇博客主要转自参考文献【1】，在原文的基础上，为了更容易增加理解，略有删改增。
一 、 重新寻找合适的细致平稳条件
在《MCMC(二)：MCMC采样和M-H采样》中，我们讲到了细致平稳条件：如果非周期马尔科夫链的状态转移矩阵 P P P和概率分布 π ( x ) \pi(x) π(x)对于所有的 i , j i,j i,j满足：
π ( i ) P ( i , j ) = π ( j ) P ( j , i ) \pi(i)P(i,j) = \pi(j)P(j,i) π(i)P(i,j)=π(j)P(j,i)
则称概率分布 π ( x ) \pi(x) π(x)是状态转移矩阵 P P P的平稳分布。
在M-H采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。
从二维的数据分布开始，假设 π ( x 1 , x 2 ) \pi(x_1,x_2) π(x1,x2)是一个二维联合数据分布，观察第一个特征维度相同的两个点 A ( x 1 ( 1 ) , x 2 ( 1 ) ) A(x_1^{(1)},x_2^{(1)}) A(x1(1),x2(1))和 B ( x 1 ( 1 ) , x 2 ( 2 ) ) B(x_1^{(1)},x_2^{(2)}) B(x1(1),x2(2))，容易发现下面两式成立：
π ( x 1 ( 1 ) , x 2 ( 1 ) ) π ( x 2 ( 2 ) ∣ x 1 ( 1 ) ) = π ( x 1 ( 1 ) ) π ( x 2 ( 1 ) ∣ x 1 ( 1 ) ) π ( x 2 ( 2 ) ∣ x 1 ( 1 ) ) (1) \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)})\qquad \text{(1)} π(x1(1),x2(1))π(x2(2)∣x1(1))=π(x1(1))π(x2(1)∣x1(1))π(x2(2)∣x1(1))(1)
π ( x 1 ( 1 ) , x 2 ( 2 ) ) π ( x 2 ( 1 ) ∣ x 1 ( 1 ) ) = π ( x 1 ( 1 ) ) π ( x 2 ( 2 ) ∣ x 1 ( 1 ) ) π ( x 2 ( 1 ) ∣ x 1 ( 1 ) ) (2) \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)})\qquad \text{(2)} π(x1(1),x2(2))π(x2(1)∣x1(1))=π(x1(1))π(x2(2)∣x1(1))π(x2(1)∣x1(1))(2)
成立的原因是： π ( x 1 ( 1 ) , x 2 ( 1 ) ) = π ( x 1 ( 1 ) ) π ( x 2 ( 1 ) ∣ x 1 ( 1 ) ) \pi(x_1^{(1)},x_2^{(1)})=\pi(x_1^{(1)})\pi(x_2^{(1)}|x_1^{(1)}) π(x1(1),x2(1))=π(x1(1))π(x2(1)∣x1(1))和 π ( x 1 ( 1 ) , x 2 ( 2 ) ) = π ( x 1 ( 1 ) ) π ( x 2 ( 2 ) ∣ x 1 ( 1 ) ) \pi(x_1^{(1)},x_2^{(2)})= \pi(x_1^{(1)}) \pi(x_2^{(2)} | x_1^{(1)}) π(x1(1),x2(2))=π(x1(1))π(x2(2)∣x1(1))。
由于式（1）（2）的右边相等，因此我们有：
π ( x 1 ( 1 ) , x 2 ( 1 ) ) π ( x 2 ( 2 ) ∣ x 1 ( 1 ) ) = π ( x 1 ( 1 ) , x 2 ( 2 ) ) π ( 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)}) π(x1(1),x2(1))π(x2(2)∣x1(1))=π(x1(1),x2(2))π(x2(1)∣x1(1))
也就是：
π ( A ) π ( x 2 ( 2 ) ∣ x 1 ( 1 ) ) = π ( B ) π ( x 2 ( 1 ) ∣ x 1 ( 1 ) ) \pi(A) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(B) \pi(x_2^{(1)} | x_1^{(1)}) π(A)π(x2(2)∣x1(1))=π(B)π(x2(1)∣x1(1))
观察上式再观察细致平稳条件的公式，我们发现在 x 1 = x 1 ( 1 ) x_1 = x_1^{(1)} x1=x1(1)这条直线上，如果用条件概率分布 π ( x 2 ∣ x 1 ( 1 ) ) \pi(x_2| x_1^{(1)}) π(x2∣x1(1))作为马尔科夫链的状态转移概率，则任意两个点之间的转移满足细致平稳条件！这真是一个开心的发现，同样的道理，在 x 2 = x 2 ( 1 ) x_2 = x_2^{(1)} x2=x2(1)这条直线上，如果用条件概率分布 π ( x 1 ∣ x 2 ( 1 ) ) \pi(x_1| x_2^{(1)}) π(x1∣x2(1))作为马尔科夫链的状态转移概率，则任意两个点之间的转移也满足细致平稳条件。那是因为假如有一点 C ( x 1 ( 2 ) , x 2 ( 1 ) ) C(x_1^{(2)},x_2^{(1)}) C(x1(2),x2(1)),我们可以得到：
π ( A ) π ( x 1 ( 2 ) ∣ x 2 ( 1 ) ) = π ( C ) π ( x 1 ( 1 ) ∣ x 2 ( 1 ) ) \pi(A) \pi(x_1^{(2)} | x_2^{(1)}) = \pi(C) \pi(x_1^{(1)} | x_2^{(1)}) π(A)π(x1(2)∣x2(1))=π(C)π(x1(1)∣x2(1))
基于上面的发现，我们可以这样构造分布 π ( x 1 , x 2 ) \pi(x_1,x_2) π(x1,x2)的马尔可夫链对应的状态转移矩阵 P P P：
P ( A → B ) = π ( x 2 ( B ) ∣ x 1 ( 1 ) )      i f    x 1 ( A ) = x 1 ( B ) = x 1 ( 1 ) (3) P(A \to B) = \pi(x_2^{(B)}|x_1^{(1)})\;\; if\; x_1^{(A)} = x_1^{(B)} =x_1^{(1)}\qquad \text{(3)} P(A→B)=π(x2(B)∣x1(1))ifx1(A)=x1(B)=x1(1)(3)
P ( A → C ) = π ( x 1 ( C ) ∣ x 2 ( 1 ) )      i f    x 2 ( A ) = x 2 ( C ) = x 2 ( 1 ) (4) P(A \to C) = \pi(x_1^{(C)}|x_2^{(1)})\;\; if\; x_2^{(A)} = x_2^{(C)} =x_2^{(1)}\qquad \text{(4)} P(A→C)=π(x1(C)∣x2(1))ifx2(A)=x2(C)=x2(1)(4)
P ( A → D ) = 0      e l s e P(A \to D) = 0\;\; else P(A→D)=0else
有了公式（3）（4）这样的状态转移矩阵P，我们很容易验证平面上的任意两点 E , F E,F E,F，满足细致平稳条件：
π ( E ) P ( E → F ) = π ( F ) P ( F → E ) \pi(E)P(E \to F) = \pi(F)P(F \to E) π(E)P(E→F)=π(F)P(F→E)
二、二维Gibbs采样
利用上一节找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样必须要知道两个维度之间的条件概率。具体过程如下：
1）输入平稳分布 π ( x 1 , x 2 ) \pi(x_1,x_2) π(x1,x2)，设定状态转移次数阈值 n 1 n_1 n1，需要的样本个数 n 2 n_2 n2
2）随机初始化初始状态值 x 1 ( 0 ) x_1^{(0)} x1(0)和 x 2 ( 0 ) x_2^{(0)} x2(0)
3）for t = 0 t=0 t=0 to n 1 + n 2 − 1 n_1+n_2−1 n1+n2−1:
a) 从条件概率分布 P ( x 2 ∣ x 1 ( t ) ) P(x_2|x_1^{(t)}) P(x2∣x1(t))中采样得到样本 x 2 t + 1 x_2^{t+1} x2t+1
b) 从条件概率分布 P ( x 1 ∣ x 2 ( t + 1 ) ) P(x_1|x_2^{(t+1)}) P(x1∣x2(t+1))中采样得到样本 x 1 t + 1 x_1^{t+1} x1t+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^{(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)})\} {(x1(n1),x2(n1)),(x1(n1+1),x2(n1+1)),...,(x1(n1+n2−1),x2(n1+n2−1))}即为我们需要的平稳分布对应的样本集。
整个采样过程中，我们通过轮换坐标轴，采样的过程为：
( x 1 ( 1 ) , x 2 ( 1 ) ) → ( x 1 ( 1 ) , x 2 ( 2 ) ) → ( x 1 ( 2 ) , x 2 ( 2 ) ) → . . . → ( 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)}) (x1(1),x2(1))→(x1(1),x2(2))→(x1(2),x2(2))→...→(x1(n1+n2−1),x2(n1+n2−1))
用下图可以很直观的看出，采样是在两个坐标轴上不停的轮换的。当然，坐标轴轮换不是必须的，我们也可以每次随机选择一个坐标轴进行采样。不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
三、多维Gibbs采样
上面的这个算法推广到多维的时候也是成立的。比如一个 n n n维的概率分布 π ( x 1 , x 2 , . . . x n ) \pi(x_1,x_2,...x_n) π(x1,x2,...xn)，我们可以通过在 n n n个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴 x i x_i xi上的转移，马尔科夫链的状态转移概率为 P ( x i ∣ x 1 , x 2 , . . . , x i − 1 , x i + 1 , . . . , x n ) P(x_i|x_1,x_2,...,x_{i-1},x_{i+1},...,x_n) P(xi∣x1,x2,...,xi−1,xi+1,...,xn)，即固定 n − 1 n−1 n−1个坐标轴，在某一个坐标轴上移动。
具体的算法过程如下:
1）输入平稳分布 π ( x 1 , x 2 ， . . . , x n ) \pi(x_1,x_2，...,x_n) π(x1,x2，...,xn)或者对应的所有特征的条件概率分布，设定状态转移次数阈值 n 1 n_1 n1，需要的样本个数 n 2 n_2 n2
2）随机初始化初始状态值 ( x 1 ( 0 ) , x 2 ( 0 ) , . . . , x n ( 0 ) ) (x_1^{(0)},x_2^{(0)},...,x_n^{(0)}) (x1(0),x2(0),...,xn(0))
3）for t = 0 t=0 t=0 to n 1 + n 2 − 1 n_1+n_2−1 n1+n2−1:
a) 从条件概率分布 P ( x 1 ∣ x 2 ( t ) , x 3 ( t ) , . . . , x n ( t ) ) P(x_1|x_2^{(t)}, x_3^{(t)},...,x_n^{(t)}) P(x1∣x2(t),x3(t),...,xn(t))中采样得到样本 x 1 t + 1 x_1^{t+1} x1t+1
b) 从条件概率分布 P ( x 2 ∣ x 1 ( t + 1 ) , x 3 ( t ) , x 4 ( t ) , . . . , x n ( t ) ) P(x_2|x_1^{(t+1)}, x_3^{(t)}, x_4^{(t)},...,x_n^{(t)}) P(x2∣x1(t+1),x3(t),x4(t),...,xn(t))中采样得到样本 x 2 t + 1 x_2^{t+1} x2t+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 ) ) P(x_j|x_1^{(t+1)}, x_2^{(t+1)},..., x_{j-1}^{(t+1)},x_{j+1}^{(t)}...,x_n^{(t)}) P(xj∣x1(t+1),x2(t+1),...,xj−1(t+1),xj+1(t)...,xn(t))中采样得到样本 x j t + 1 x_j^{t+1} xjt+1
f)…
g) 从条件概率分布 P ( x n ∣ x 1 ( t + 1 ) , x 2 ( t + 1 ) , . . . , x n − 1 ( t + 1 ) ) P(x_n|x_1^{(t+1)}, x_2^{(t+1)},...,x_{n-1}^{(t+1)}) P(xn∣x1(t+1),x2(t+1),...,xn−1(t+1))中采样得到样本 x n t + 1 x_n^{t+1} xnt+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 ) ) } \{(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)})\} {(x1(n1),x2(n1),...,xn(n1)),...,(x1(n1+n2−1),x2(n1+n2−1),...,xn(n1+n2−1))}即为我们需要的平稳分布对应的样本集。
整个采样过程和Lasso回归的坐标轴下降法算法非常类似，只不过Lasso回归是固定 n − 1 n−1 n−1个特征，对某一个特征求极值。而Gibbs采样是固定 n − 1 n−1 n−1个特征在某一个特征采样。
同样的，轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
四、二维Gibbs采样实例
这里给出一个Gibbs采样的例子。
假设我们要采样的是一个二维正态分布 N o r m ( μ , Σ ) Norm(\mu,\Sigma) Norm(μ,Σ),其中：
μ = ( μ 1 , μ 2 ) = ( 5 , − 1 ) \mu = (\mu_1,\mu_2) = (5,-1) μ=(μ1,μ2)=(5,−1)
Σ = ( σ 1 2 ρ σ 1 σ 2 ρ σ 1 σ 2 σ 2 2 ) = ( 1 1 1 4 ) \Sigma = \left( \begin{array}{ccc} \sigma_1^2&\rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 &\sigma_2^2 \end{array} \right) = \left( \begin{array}{ccc} 1&1 \\ 1&4 \end{array} \right) Σ=(σ12ρσ1σ2ρσ1σ2σ22)=(1114)
而采样过程中的需要的状态转移条件分布即为二维正态分布的条件概率分布为：
P ( x 1 ∣ x 2 ) = N o r m ( μ 1 + ρ σ 1 / σ 2 ( x 2 − μ 2 ) , ( 1 − ρ 2 ) σ 1 2 ) P(x_1|x_2) = Norm\left ( \mu _1+\rho \sigma_1/\sigma_2 \left ( x _2-\mu _2 \right ), (1-\rho ^2)\sigma_1^2 \right ) P(x1∣x2)=Norm(μ1+ρσ1/σ2(x2−μ2),(1−ρ2)σ12)
P ( x 2 ∣ x 1 ) = N o r m ( μ 2 + ρ σ 2 / σ 1 ( x 1 − μ 1 ) , ( 1 − ρ 2 ) σ 2 2 ) P(x_2|x_1) = Norm\left ( \mu _2+\rho \sigma_2/\sigma_1 \left ( x _1-\mu _1 \right ), (1-\rho ^2)\sigma_2^2 \right ) P(x2∣x1)=Norm(μ2+ρσ2/σ1(x1−μ1),(1−ρ2)σ22)
具体的代码如下：
import random
import math
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
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 range(N):
for j in range(K):
x = p_xgiveny(y, m1, m2, s1, s2)
y = p_ygivenx(x, m1, m2, s1, s2)
z = samplesource.pdf([x,y])
# x_res,y_res,z_res的长度均为100,000.
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()
输出的正态分布图如下：
这里对代码说明一点：
①代码里双重for循环的用法的确和算法描述稍有不同：
在代码里外层循环是得到 Ｎ Ｎ Ｎ个样本，对于每个样本，我们要在内层循环中分别进行马尔科夫状态转移 Ｋ Ｋ Ｋ次得到对应的样本。这样的好处是得到的样本独立性好一点，毕竟我们的条件概率计算量不大。
在算法里描述得是，先进行若干次马尔科夫状态转移，到了我们认为已经收敛到目标分布时，就进行连续的采样，得到 Ｎ Ｎ Ｎ个样本。这样的好处是采样的计算量小很多，尤其是对于复杂分布的计算，一次马尔科夫过程后就可以得到所有想要个数的样本，但是样本之间的独立性会稍差。
五、Gibbs采样小结
由前几节我们知道：Gibbs采样必须要知道特征维度之间的条件概率。如果不知道这个转移的条件概率，则无法进行Gibbs采样。
同时由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。
有了Gibbs采样来获取概率分布的样本集，有了蒙特卡罗方法来用样本集模拟求和，他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。
这里，我们对3篇MCMC的博客进行一个总的结论，3篇博客的最大作用是：我们得到了正确采样的方法论，当有了这个方法后，一些复杂的问题比如基于概率分布的积分求和问题可以使用采样的方法近似求解。
比如在LDA主题模型中，我们可以用MCMC来采样求主题分布和词分布的近似解。
文本主题模型之LDA(二) LDA求解之Gibbs采样算法
此外还有分解机模型里的积分求解。
分解机(Factorization Machines)推荐算法原理
参考文献
【1】MCMC(四)Gibbs采样