gibbs sampling

  原文链接：http://www.algorithmdog.com/mh-and-gibbs?utm_source=tuicool&utm_medium=referral

在科学研究中，如何生成服从某个概率分布的样本是一个重要的问题。 如果样本维度很低，只有一两维，我们可以用反切法、拒绝采样和重要性采样等方法。 但是对于高维样本，这些方法就不适用了。这时我们就要使用一些 “高档” 的算法，比如下面要介绍的 Metropolis-Hasting 算法和 Gibbs sampling 算法。

      Metropolis-Hasting 算法和 Gibbs sampling 算法是马尔科夫链蒙特卡洛（Markov Chain Mento Carlo, MCMC）方法。我们先介绍 MCMC 方法。

1. 马尔科夫蒙特卡洛方法

      MCMC 方法是用蒙特卡洛方法去体现马尔科夫链的方法。马尔科夫链是状态空间的转换关系，下一个状态只和当前的状态有关。比如下图就是一个马尔科夫链的示意图。

图中转移关系可以用一个概率转换矩阵 p 表示，

p=⎡⎣⎢000.610.10.400.90⎤⎦⎥(1)

如果当前状态分布为 u(x), 那么下一个矩阵的状态就是 u(x)p, 再下一个就是u(x)p2,... 最后会收敛到一个平稳分布 π。这个平稳分布 π 只和概率转移矩阵 p 有关，而和初始状态分布 u 是什么没有关系。

      如何判断一个马尔科夫链是否能收敛到平稳分布，以及如何判断一个状态分布是不是一个马尔科夫链的平稳分布呢？我们有下面定理。

细致平衡条件: 已知各态历经的的马尔科夫链有概率转移矩阵 p，以及已知状态分布 π。如果对于任意两个状态 i 和 j，下面公式成立，则马尔科夫链能够收敛到 π。

π(i)p(j|i)=π(j)p(i|j)

这里的各态历经是指任意两个状态之间可以通过有限步到达。

怎么证明细致平衡条件呢？我也不知道啊。

      MCMC 方法的基本原理是利用细致平衡条件构建一个概率转移矩阵，使得目标概率就是概率转移矩阵的平稳分布。 Metropolis-Hasting 和 Gibbs sampling 算法本质上是构建概率转移矩阵的不同方法。

2. Metropolis-Hastings 算法

      Metropolis-Hastings 算法先提出一个可能不符合条件的概率转移矩阵 q, 然后再进行调整。比如我们提出的 q 是均匀概率，即从任意状态到任意状态的概率是相等的。显然在绝大部分情况下，q 的稳定概率不是目标概率 π，即不满足细致平衡条件。

π(i)q(j|i)≠π(j)q(i|j)

如何让这个不等式转变成等式呢？根据对称性，我们容易得到下面的等式。

π(i)q(j|i)π(j)q(i|j)=π(j)q(i|j)π(i)q(j|i)(2)

这时整个概率转移矩阵满足细致平衡条件。从 i 状态转到 j 状态的概率是 q(j|i)π(j)q(i|j)，实现这个转移概率的方式是 i 状态以 q(j|i) 概率跳转到 j 状态，然后以 π(j)q(i|j) 接受跳转 (拒绝跳转就退回 i 状态)。这样整个 Metropolis-Hasting 算法的框架就建立起来了。

      这个原始的 Metropoli-Hasting 算法的有一个小问题。 跳转接受概率 π(j)q(i|j) 和 π(i)q(j|i) 的值很小，算法进行过程充斥着跳转拒绝。为了改进这点，Metropoli-Hasting 算法的方法是公式两边同时乘以一个系数，使得 π(j)q(i|j) 和 π(i)q(j|i) 中大的一项 scale 到 1，得到下面的公式。

π(i)q(j|i)π(j)q(i|j)π(i)q(j|i)orπ(i)q(j|i)==π(j)q(i|j)whenπ(i)q(j|i)>π(j)q(i|j)π(j)q(i|j)π(i)q(j|i)π(j)q(i|j)whenπ(i)q(j|i)≤π(j)q(i|j)

这个公式可以进一步简化为下面的公式

π(i)q(j|i)a(j|i)a(j|i)a(i|j)===π(j)q(i|j)a(i|j)min{π(j)q(i|j)π(i)q(j|i),1}min{π(i)q(j|i)π(j)q(i|j),1}(3)

      根据上面的推导，我们容易得到 Metropolis-Hasting 算法的流程。

3. Gibbs sampling 算法

      Gibbs sampling 算法是 Metropolis-Hasting 算法的一个特例。很鸡贼的一个特例。m 维的一个样本跳转到另一个样本的过程，可以拆解为 m 个子过程，每一个子过程对应一个维度。这时概率转移矩阵是 m 个子概率转移矩阵之积，即 p=∏mi=kpk

其中 pk 表示第 k 维的变化概率。在 pk 中，两个状态之间只有 k 维不同，其跳转概率如下所示；不然为 0。

pk(xx⊣k,k=v2|xx⊣k,k=v1)=π(xx⊣k,k=v2)∑vπ(xx⊣k,k=v)

其中 xx⊣k,k=v2 表示样本第 k 维数据为 v2，其它维度固定。这时候我们发现如下公式

===π(xx⊣k,k=v1)p(xx⊣k,k=v2|xx⊣k,k=v1)π(xx⊣k,k=v1)π(xx⊣k,k=v2)∑vπ(xx⊣k,k=v)π(xx⊣k,k=v2)π(xx⊣k,k=v1)∑vπ(xx⊣k,k=v)π(xx⊣k,k=v2)p(xx⊣k,k=v1|xx⊣k,k=v2)

即 pk 和 π 满足细致平衡条件的等式。那么 pk 就是我们要构建的概率转移矩阵嘛？答案是否定的。因为完整的细致平衡条件需要各态历经。在概率转移矩阵 pk 下, 只有 k 维数据子啊变化，因此一个状态永远不能到达和它第 k-1 维数据不同的状态。

      最终我们构建的概率转移矩阵是 m 个子概率转移矩阵之积

p=∏i=kmpk(4)

我们很容易证明 p 依然满足细致平衡条件中的等式，同时还满足各态历经。根据上述推导，我们得到 Gibbs sampling 的算法过程。

4. 总结

      Metropolist-Hasting 和 Gibbs sampling 是有效的 MCMC 算法，能够解决高维空间的采样问题。