转载

Metropolis-Hastings 算法和 Gibbs sampling 算法 | AlgorithmDog

我们有一个概率分布 /(/pi/)

$/pi$ ，需要生成满足这个分布的样本。如果样本维度很低，只有一两维，我们可以用反切法、拒绝采样和重要性采样等方法。但是样本维度很高，成千上百，这些方法就不适用了。这时我们就要使用一些高档的算法，比如下面要介绍的 Metropolis-Hasting 算法和 Gibbs sampling 算法。

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

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

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

Metropolis-Hastings 算法和 Gibbs sampling 算法 | AlgorithmDog

图中转移关系可以用一个概率转换矩阵 p 表示，
/begin{eqnarray}
P = /begin{bmatrix}
0 &1 &0 //
0&0.1 &0.9 //
0.6 &0.4 &0
/end{bmatrix}
/end{eqnarray}

如果当前状态分布为 /(u(x) = (0.3,0.2,0.5)/)

$u(x) = (0.3,0.2,0.5)$ , 那么下一个矩阵的状态就是 /( u(x)p /) $u(x)p$ , 再下一个就是 /(u(x)p^2/) $u(x)p^2$ ,... 最后会收敛到一个平稳分布 /(/pi/) $/pi$ 。这个平稳分布 /(/pi/) $/pi$ 只和概率转移矩阵 p 有关，而和初始状态分布 u 是什么没有关系。

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

细致平衡条件: 已知各态历经的的马尔科夫链有概率转移矩阵 /(p/)

$p$ ，以及已知状态分布 /(/pi/) $/pi$ 。如果对于任意两个状态 i 和 j，下面公式成立，则马尔科夫链能够收敛到 /(/pi/) $/pi$ 。
/begin{eqnarray}
/pi(i) p(j|i) = /pi(j) p(i|j) /nonumber
/end{eqnarray}
这里的各态历经是指任意两个状态之间可以通过有限步到达。

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

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

2. Metropolis-Hastings 算法

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

$/pi$ ，即不满足细致平衡条件。
/begin{eqnarray}
/pi(i) q(j|i) /ne /pi(j)q(i|j) /nonumber
/end{eqnarray}

如何让这个不等式转变成等式呢？根据对称性，我们容易得到下面的等式。
/begin{eqnarray}
/label{mh}
/pi(i) q(j|i) /pi(j)q(i|j) = /pi(j)q(i|j) /pi(i)q(j|i)
/end{eqnarray}

这时整个概率转移矩阵满足细致平衡条件。从 i 状态转到 j 状态的概率是 /(q(j|i) /pi(j)q(i|j)/)

这个原始的 Metropoli-Hasting 算法的有一个小问题。跳转接受概率 /(/pi(j)q(i|j)/)

/begin{eqnarray}
/pi(i) q(j|i) /frac{/pi(j)q(i|j)}{/pi(i)q(j|i)} &=& /pi(j)q(i|j) /;/; when /;/; /pi(i)q(j|i) > /pi(j)q(i|j) /nonumber //
or/nonumber //
/pi(i) q(j|i) &=& /pi(j)q(i|j)/frac{/pi(i)q(j|i)}{/pi(j)q(i|j)} /;/; when /;/; /pi(i)q(j|i) /le /pi(j)q(i|j) /nonumber //
/end{eqnarray}

这个公式可以进一步简化为公式 2
/begin{eqnarray}
/pi(i) q(j|i) a(j|i) &=& /pi(j)q(i|j) a(i|j) /nonumber //
a(j|i) &=& min/{/frac{/pi(j)q(i|j)}{/pi(i)q(j|i)},1/} /nonumber //
a(i|j) &=& min/{/frac{/pi(i)q(j|i)}{/pi(j)q(i|j)},1/}
/end{eqnarray}

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

3. Gibbs sampling 算法

Gibbs sampling 算法是 Metropolis-Hasting 算法的一个特例。很鸡贼的一个特例。m 维的一个样本跳转到另一个样本的过程，可以拆解为 m 个子过程，每一个子过程对应一个维度。这时概率转移矩阵可以看成是 m 个子概率转移矩阵之积，即 /(p = /prod_{i=k}^{m} p_k /)

$p = /prod_{i=k}^{m} p_k$
/end{eqnarray}

其中 /(p_k/)

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

/begin{eqnarray}
p_k(/pmb{x}_{/dashv k, k=v_2}|/pmb{x}_{/dashv k, k=v_1}) = /frac{/pi(/pmb{x}_{/dashv k, k=v_2})}{/sum_{v}/pi(/pmb{x}_{/dashv k, k=v})} /nonumber
/end{eqnarray}

其中 /(/pmb{x}_{/dashv k, k=v_2}/)

$/pmb{x}_{/dashv k, k=v_2}$ 表示样本第 k 维数据为 /(v_2/) $v_2$ ，其他没有变化。这时候我们发现如下公式

/begin{eqnarray}
&& /pi(/pmb{x}_{/dashv k, k=v_1}) p(/pmb{x}_{/dashv k, k=v_2}|/pmb{x}_{/dashv k, k=v_1}) /nonumber //
&=& /pi(/pmb{x}_{/dashv k, k=v_1}) /frac{/pi(/pmb{x}_{/dashv k, k=v_2})}{/sum_{v}/pi(/pmb{x}_{/dashv k, k=v})} /nonumber //
&=& /pi(/pmb{x}_{/dashv k, k=v_2}) /frac{/pi(/pmb{x}_{/dashv k, k=v_1})}{/sum_{v}/pi(/pmb{x}_{/dashv k, k=v})} /nonumber //
&=& /pi(/pmb{x}_{/dashv k, k=v_2}) p(/pmb{x}_{/dashv k, k=v_1}|/pmb{x}_{/dashv k, k=v_2}) /nonumber //
/end{eqnarray}
即 /(p_k/)

$p_k$ 和 /(/pi/) $/pi$ 满足细致平衡条件的等式。那么 /(p_k/) $p_k$ 就是我们要构建的概率转移矩阵嘛？答案是否定的。因为完整的细致平衡条件需要各态历经，而一个状态在概率转移矩阵 /(p_k/) $p_k$ 下永远不能到达另一个第 k-1 维数据不同的状态。我们最终构建的概率转移矩阵
/begin{eqnarray}
p = /prod_{i=k}^{m} p_k
/end{eqnarray}

我们很容易证明 /(p/)

$p$ 依然满足细致平衡条件中的等式，同时还满足各态历经。根据这些推理，我们得到 Gibbs sampling 的算法过程。