引言

上一节介绍了马尔可夫链(Markov Chain)以及平稳分布。本节将马尔可夫链与蒙特卡洛方法相结合，介绍MH采样算法。

回顾：马尔可夫链与平稳分布

马尔可夫链

基于齐次马尔可夫假设，以一阶齐次马尔可夫假设 为例，$t+1$时刻的随机变量${\mathcal{X}}_{t+1}$与$t$时刻的随机变量${\mathcal{X}}_t$之间的关联关系如下：
$$P({\mathcal{X}}_{t+1}\mid {\mathcal{X}}_t,{\mathcal{X}}_{t-1},\cdots ,{\mathcal{X}}_1)=P({\mathcal{X}}_{t+1}\mid {\mathcal{X}}_t)$$
针对马尔可夫链中任一时刻的随机变量${\mathcal{X}}_t$可选择的结果均是离散的 这种情况，定义共包含$\mathcal{K}$种选择方式，$t$时刻随机变量${\mathcal{X}}_t$的概率分布$\pi ({\mathcal{X}}_t)$表示如下：
$$\pi ({\mathcal{X}}_t)={\left[\pi ({\mathcal{X}}_t=x_1),\pi ({\mathcal{X}}_t=x_2),\cdots ,\pi ({\mathcal{X}}_t=x_{\mathcal{K}})\right]}_{\mathcal{K}\times 1}^T$$
对应地，该马尔可夫链的状态转移矩阵是一个$\mathcal{K}\times \mathcal{K}$的方阵；并且矩阵中的每一个元素$a_{ij}\in \mathcal{A}$均表示基于转移过程选择的条件概率：
A = [ a i j ] K × K i , j ∈ { 1 , 2 , ⋯   , K } a i j = P ( X t + 1 = x j ∣ X t = x i ) t = 1 , 2 , ⋯ \mathcal A = [a_{ij}]_{\mathcal K \times \mathcal K} \quad i,j \in \{1,2,\cdots,\mathcal K\} \\ a_{ij} = P(\mathcal X_{t+1} = x_j \mid \mathcal X_t = x_i) \quad t=1,2,\cdots A=[aij​]K×K​i,j∈{1,2,⋯,K}aij​=P(Xt+1​=xj​∣Xt​=xi​)t=1,2,⋯

平稳分布

平稳分布(Stationary Distribution)表示马尔可夫链具有某种平稳性质的概率分布。平稳分布存在的底层逻辑在于 状态转移矩阵是一个双随机矩阵：

• 状态转移矩阵是一个$\mathcal{K}\times \mathcal{K}$的方阵；
• 状态转移矩阵各行、列元素之和为1；

基于双随机矩阵的性质，得到状态转移矩阵的特征值均$\le 1$恒成立。
通过证明，只要马尔可夫链的状态转移矩阵 确定的条件下，在未来的转移过程中必然会逼近至平稳分布。
具体证明过程请看上一节内容~传送门

假设$m$时刻状态下达到平稳分布，则有：
$$\pi ({\mathcal{X}}_m)=\pi ({\mathcal{X}}_{m+1})=\pi ({\mathcal{X}}_{m+2})=\cdots$$

如何判定当前时刻的分布是否为平稳分布？需要借助细致平衡原则(Detail Balance)：
$$\pi (\mathcal{X}=x_i)\cdot P(\mathcal{X}=x^{\ast }\mid \mathcal{X}=x)=\pi (\mathcal{X}=x^{\ast })\cdot P(\mathcal{X}=x\mid \mathcal{X}=x^{\ast })$$
细致平衡是概率分布是平稳分布的充分非必要条件，因此可以通过细致平衡去判别平稳分布，反之不然。

MH采样算法

在蒙特卡洛方法介绍中提到，推断(Inference)关心的问题是后验概率$P(\mathcal{Z}\mid \mathcal{X})$的结果，或者是关于$P(\mathcal{Z}\mid \mathcal{X})$的期望信息：
$${\mathbb{E}}_{\mathcal{Z}\mid \mathcal{X}}[f(\mathcal{Z})]={\int }_{\mathcal{Z}\mid \mathcal{X}}P(\mathcal{Z}\mid \mathcal{X})f(\mathcal{Z})d\mathcal{Z}$$
通过蒙特卡洛方法，从概率分布$P(\mathcal{Z}\mid \mathcal{X})$中进行采样，再进行计算：
$$\begin{gathered} z^{(1)},z^{(2)},\cdots ,z^{(N)}\sim P(\mathcal{Z}\mid \mathcal{X}) \\ {\mathbb{E}}_{\mathcal{Z}\mid \mathcal{X}}[f(\mathcal{Z})]=\frac{1}{N}\sum _{i=1}^{N}f(z^{(i)}) \end{gathered}$$
但真实情况是：从$P(\mathcal{Z}\mid \mathcal{X})$中采集样本是非常复杂的事情。因此借助马尔可夫链(Markov Chain)来间接获取样本信息。

采样思路

马尔可夫链是如何实现采样的？
假设已经得到一个状态转移矩阵${\mathcal{A}}^{\ast }$，${\mathcal{A}}^{\ast }$满足：在一阶齐次马尔可夫假设的条件下，任意两个连续状态下的随机变量$\mathcal{X}$对应的概率分布，均满足细致平衡原则。

换句话说，${\mathcal{A}}^{\ast }$使马尔可夫链 { X T } \{\mathcal X_{T}\} {XT​}共用同一个概率分布，也就是平稳分布。

• 此时，给定一个初始节点：
  $${\mathcal{X}}_{init}=x_i(i=1,2,\cdots ,\mathcal{K})$$
• 通过状态转移过程，随机得到下一状态随机变量的选择结果$x_j$：
  $$x_j\sim P({\mathcal{X}}_1\mid {\mathcal{X}}_{init}=x_i)$$
• 同上，根据上一时刻的状态转移过程，随机得到下一状态随机变量的选择结果$x_k$：
  $$x_k\sim P({\mathcal{X}}_2\mid {\mathcal{X}}_1=x_j)$$

以此类推，最终会得到这样一组样本点集合。这些样本点均服从平稳分布：
{ x i , x j , x k , ⋯   } \{x_i,x_j,x_k,\cdots\} {xi​,xj​,xk​,⋯}

MH采样算法过程

基于上述采样思路，将获取平稳分布问题转化为：如何找到一个恰当的状态转移矩阵，使得马尔可夫链的各分布是平稳分布。

• 首先，构建关于$P(\mathcal{Z}\mid \mathcal{X})$的马尔可夫链，并随机构建一个状态转移矩阵$\mathcal{Q}$：
  Q = [ q i j ] K × K ( i , j ∈ { 1 , 2 , ⋯   , K } ) \mathcal Q = [q_{ij}]_{\mathcal K \times \mathcal K} \quad (i,j \in \{1,2,\cdots,\mathcal K\}) Q=[qij​]K×K​(i,j∈{1,2,⋯,K})
  此时关于细致平衡的等式两项有：
  不能说是一定不相等，而是相等的概率极低。任意随机一个$\mathcal{Q}$就能得到平稳分布，那运气可太好了~
  $$P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)\ne P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })$$

• 基于上述式子，假设存在关于$z,z^{\ast }$的函数$\alpha (z,z^{\ast })$，使得：
  $$P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)\cdot \alpha (z,z^{\ast })=P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })\cdot \alpha (z^{\ast },z)$$
  我们将$Q(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)\cdot \alpha (z,z^{\ast })$记作$\mathcal{P}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)$；
  反之，$\mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })\cdot \alpha (z^{\ast },z)$记作$\mathcal{P}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })$。上述公式可转化为：
  只是一个符号的变换，$\mathcal{P}$和$P$不是同一个符号~
  $$P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{P}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)=P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{P}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })$$
  此时，上述公式满足细致平衡原则，此时的分布是平稳分布。

• 回溯上述过程，我们称$\alpha (z,z^{\ast })$函数为接收率，其具体表示如下：
  $$\alpha (z,z^{\ast })=\min \left[1,\frac{P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })}{P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)}\right]$$
  对接受率函数$\alpha (z,z^{\ast })$进行分析：

  • 观察$\alpha (z,z^{\ast })$，因为$\mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })$和$\mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })$均是条件概率，因此分数项大于等于0恒成立。从而$\alpha (z,z^{\ast })$的值域为：
    $$\alpha (z,z^{\ast })\in [0,1]$$
  • 将$\alpha (z,z^{\ast })$代入上式中，有：
    $$\begin{array}{rl}& P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)\cdot \alpha (z,z^{\ast })\\ & =P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)\cdot \min \left[1,\frac{P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })}{P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)}\right]\\ & =\min \left[P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z),P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })\right]\\ & =P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })\cdot \min \left[1,\frac{P(\mathcal{Z}=z\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)}{P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })}\right]\\ & =P(\mathcal{Z}=z^{\ast }\mid \mathcal{X})\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })\cdot \alpha (z^{\ast },z)\end{array}$$

该接收率思路即MH采样算法(Metropolis Hastings)。
它的具体算法表示如下：

• 从(0,1)均匀分布中进行采样；
  $$u\sim \mathcal{U}(0,1)$$
• 以上一时刻采样结果$z^{(t-1)}$确定的条件下，对状态转移矩阵当前时刻$z^{(t)}$进行采样；
  此时，随机采出一个结果——$z^{\ast }$,但不能直接作为$z^{(t)}$，需要满足后续的判断条件。
  $$z^{\ast }\sim \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{(t-1)})$$
• 计算$\alpha$函数的具体结果：
  $$\alpha =\min \left[1,\frac{P(\mathcal{Z}=z^{\ast })\cdot \mathcal{Q}(\mathcal{Z}=z\mid \mathcal{Z}=z^{\ast })}{P(\mathcal{Z}=z)\cdot \mathcal{Q}(\mathcal{Z}=z^{\ast }\mid \mathcal{Z}=z)}\right]$$
• 对$u$和$\alpha$结果进行如下对比：
  类似于‘拒绝采样’，从0-1均匀分布中采样的$u$自身没有任何实际意义，它只是用来衡量$\alpha$结果的性能。
  • 如果：$u\le \alpha \to z^{(t)}=z^{\ast }$；
    存在$\alpha$的概率接收$z^{\ast }$样本。
  • 否则：$z^{(t)}=z^{(t-1)}$
    此次的采样被拒绝，依然使用上一时刻的样本结果$z^{(t-1)}$作为本时刻的输出样本。这次迭代确实白跑~样本数量没有减少，下次迭代依然将$z^{(t-1)}$作为条件，基于该条件的概率分布进行采样。
• 最终会得到一系列样本：
  { z ( 1 ) , z ( 2 ) , ⋯   , z ( N ) } \{z^{(1)},z^{(2)},\cdots,z^{(N)}\} {z(1),z(2),⋯,z(N)}
  并最后通过上述样本点使用蒙特卡洛方法近似求解概率分布$P(\mathcal{Z}\mid \mathcal{X})$的期望结果：
  $${\mathbb{E}}_{\mathcal{Z}\mid \mathcal{X}}[f(\mathcal{Z})]=\frac{1}{N}\sum _{i=1}^{N}f(z^{(i)})$$

下一节将介绍吉布斯采样算法(Gibbs)

相关参考：
机器学习-白板推导系列-蒙特卡洛方法3（Monte Carlo Method）-MH采样（Metropolis Hastings）