Markov Chain Monte Carlo 采样算法

本文最后更新于：2024年5月7日下午

作为一种随机采样方法，马尔科夫链蒙特卡罗（Markov Chain Monte Carlo，以下简称MCMC）在机器学习,深度学习以及自然语言处理等领域都有广泛的应用，是很多复杂算法求解的基础，本文介绍基本思想。

简介

马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo)，简称MCMC，产生于20世纪50年代早期，是在贝叶斯理论框架下，通过计算机进行模拟的蒙特卡洛方法(Monte Carlo)。该方法将马尔科夫(Markov)过程引入到Monte Carlo模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡罗积分只能静态模拟的缺陷。MCMC是一种简单有效的计算方法，在很多领域到广泛的应用，如统计物、贝叶斯(Bayes)问题、计算机问题等。

——百度百科

背景

给定连续随机变量 $ X \in \mathcal{X} $ 的概率密度函数 $ \tilde{p}(x) $, 则 $ X $ 在区间 $ \mathbb{A} $ 中的概率可以计算为：
$$
P(\mathbb{A})=\int_{\mathrm{A}} \tilde{p}(x) d x
$$
如果函数 $ f: \mathcal{X} \longmapsto \mathbb{R} $, 则可以计算 $ f(X) $ 的期望:

$$
\mathbb{E}_{X \sim \tilde{p}(x)}[f(X)]=\int f(x) \tilde{p}(x) d x
$$

如果 $ X $ 不是一个单变量, 页是一个高维的多元变量 $ \vec{X} $, 且服从一个非常复杂的分布, 则对于上式的求积分非常困难。为此，MCMC 先构造出服从 $ \tilde{p} $ 分布的独立同分布随机变量 $ \overrightarrow{\mathbf{x}}_{1}, \overrightarrow{\mathbf{x}}_{2}, \cdots, \overrightarrow{\mathbf{x}}_{N}$ ，再得到$ \mathbb{E}_{\vec{X} \sim \tilde{p}(\overrightarrow{\mathbf{x}})}[f(\vec{X})] $ 的无偏估计:

$$ \mathbb{E}_{\vec{X} \sim \hat{p}(\vec{x})}[f(\vec{X})]=\frac{1}{N} \sum_{i=1}^{N} f\left(\overrightarrow{\mathbf{x}}_{i}\right) $$

如果概率密度函数 $ \tilde{p}(\overrightarrow{\mathbf{x}}) $ 也很复杂, 则构造服从 $ \tilde{p} $ 分布的独立同分布随机变量也很困难。 MCMC 方法就是通过构造平稳分布为 $ \tilde{p}(\overrightarrow{\mathbf{x}}) $ 的马尔可夫链来产生样本。

Metropolis Hastings 算法

基本思想

先设法构造一条马尔可夫链, 使其收敛到平稳分布恰好为 $ \tilde{p} $ 。然后通过这条马尔可夫链来产生符合 $ \tilde{p} $ 分布的样本。最后通过这些样本来进行估计。
这里马尔可夫链的构造至关重要，不同的构造方法将产生不同的 MCMC 算法。
Metropolis-Hastings : MH 算法是 MCMC 的重要代表。

构造平稳分布转移矩阵

假设已经提供了一条马尔可夫链，其转移矩阵为$Q$。目标是另一个马尔科夫链，使转移矩阵为$P$、平稳分布是 $ \tilde{p} $ 。
通常 $ \tilde{p}(i) Q_{i, j} \neq \tilde{p}(j) Q_{j, i} $, 即 $ \tilde{p} $ 并不满足细致平稳条件
我们需要改造已有的马尔可夫链，使得细致平稳条件成立
引入一个函数 $ \alpha(i, j) $, 使其满足:

$$
\tilde{p}(i) Q_{i, j} \alpha(i, j)=\tilde{p}(j) Q_{j, i} \alpha(j, i)
$$

可以取 $ \alpha(i, j)=\tilde{p}(j) Q_{j, i} $, 则有:

$$
\tilde{p}(i) Q_{i, j} \alpha(i, j)=\tilde{p}(i) Q_{i, j} \tilde{p}(j) Q_{j, i}=\tilde{p}(j) Q_{j, i} \tilde{p}(i) Q_{i, j}=\tilde{p}(j) Q_{j, i} \alpha(j, i)
$$

令： $ Q_{i, j}^{\prime}=Q_{i, j} \alpha(i, j), Q_{j, i}^{\prime}=Q_{j, i} \alpha(j, i) $, 则有 $ \tilde{p}(i) Q_{i, j}^{\prime}=\tilde{p}(j) Q_{j, i}^{\prime} $其中 $ Q_{i, j}^{\prime} $ 构成了转移矩阵 $ \mathbf{Q}^{\prime} $。而 $ \mathbf{Q}^{\prime} $ 满足细致平稳条件，因此它对应的马尔可夫链的平稳分布就是 $ \tilde{p} $ 。

接受率

在改造 $ \mathbf{Q} $ 的过程中引入的 $ \alpha(i, j) $ 称作接受率。

其物理意义为: 在原来的马尔可夫链上，从状态 $ i $ 以 $ Q_{i, j} $ 的概率跳转到状态$j$的时候，以$ \alpha(i, j) $的概率接受这个转移。
如果接受率$ \alpha(i, j) $太小，则改造马尔可夫链过程中非常容易原地踏步，拒绝大量的跳转。
这样使得马尔可夫链遍历所有的状态空间需要花费太长的时间，收敛到平稳分布$ \tilde{p} $的速度太慢。
根据推导 $ \alpha(i, j)=\tilde{p}(j) Q_{j, i} $, 如果将系数从 1 提高到 $ K $, 则有:

$$ \begin{aligned} \alpha^{*}(i, j)=K \tilde{p}(j) Q_{j, i} &=K \alpha(i, j) \\ \alpha^{*}(j, i)=K \tilde{p}(i) Q_{i, j} &=K \alpha(j, i) \end{aligned} $$

于是有:

$$ \tilde{p}(i) Q_{i, j} \alpha^{*}(i, j)=K \tilde{p}(i) Q_{i, j} \alpha(i, j)=K \tilde{p}(j) Q_{j, i} \alpha(j, i)=\tilde{p}(j) Q_{j, i} \alpha^{*}(j, i) $$

因此，即使提高了接受率, 细致平稳条件仍然成立。

将 $ \alpha(i, j), \alpha(j, i) $ 同比例放大, 取: $ \alpha(i, j)=\min \left\{\frac{\tilde{p}(j) Q_{j, i}}{\tilde{p}(i) Q_{i, j}}, 1\right\} $ 。

当$ \tilde{p}(j) Q_{j, i}=\tilde{p}(i) Q_{i, j} $ 时： $ \alpha(i, j)=\alpha(j, i)=1 $, 此时满足细致平稳条件。
当$ \tilde{p}(j) Q_{j, i}>\tilde{p}(i) Q_{i, j} $ 时： $ \alpha(i, j)=1, \alpha(j, i)=\frac{\tilde{p}(i) Q_{i, j}}{\tilde{p}(j) Q_{j, i}} $, 此时满足细致平稳条件。
当$ \tilde{p}(j) Q_{j, i}<\tilde{p}(i) Q_{i, j} $ 时： $ \alpha(i, j)=\frac{\tilde{p}(j) Q_{j, j}}{\tilde{p}(i) Q_{i, j}}, \alpha(j, i)=1 $, 此时满足细致平稳条件。

算法

输入:

先验转移概率矩阵$Q$
目标分布 $ \tilde{p} $

输出: 采样出的一个状态序列 $ \left\{x_{0}, x_{1}, \cdots, x_{n}, x_{n+1}, \cdots\right\} $

算法步骤：

初始化 $ x_{0} $
对 $ t=1,2, \cdots $ 执行迭代
迭代步骤:
- 根据$ Q\left(x^{*} \mid x_{t-1}\right) $采样出候选样本$ x^{*} $, 其中 $ Q $ 为转移概率函数。
- 计算 $ \alpha\left(x^{*} \mid x_{t-1}\right) $ :
$$ \alpha\left(x^{*} \mid x_{t-1}\right)=\min \left(1, \frac{\tilde{p}\left(x^{*}\right) Q\left(x_{t-1} \mid x^{*}\right)}{\tilde{p}\left(x_{t-1}\right) Q\left(x^{*} \mid x_{t-1}\right)}\right) $$
- 根据均匀分布从 $ (0,1) $ 内采样出阈值 $ u $, 如果 $ u \leq \alpha\left(x^{*} \mid x_{t-1}\right) $, 则接受 $ x^{*} $, 即 $ x_{t}=x^{*} $ 。否
- 返回采样得到的序列 $ \left\{x_{0}, x_{1}, \cdots, x_{n}, x_{n+1}, \cdots\right\} $
注意：返回的序列中，只有充分大的 $ n $ 之后的序列 $ \left\{x_{n}, x_{n+1}, \cdots\right\} $ sss才是服从 $ \tilde{p} $ 分布的采样序列。

Gibbs 算法

MH 算法不仅可以应用于一维空间的采样，也适合高维空间的采样。
对于高维的情况，由于接受率 $ \alpha $ 的存在 (通常 $ \alpha<1) $ ， MH 算法的效率通常不够高, 此时一般使用 Gibbs 采样算法。

考虑二维的情形：假设有概率分布$ \tilde{p}(x, y) $, 考察状态空间上 $ x $ 坐标相同的两个点 $ A\left(x_{1}, y_{1}\right), B\left(x_{1}, y_{2}\right) $，可以证明有:

$$ \begin{array}{l} \tilde{p}\left(x_{1}, y_{1}\right) \tilde{p}\left(y_{2} \mid x_{1}\right)=\tilde{p}\left(x_{1}\right) \tilde{p}\left(y_{1} \mid x_{1}\right) \tilde{p}\left(y_{2} \mid x_{1}\right) \\ \tilde{p}\left(x_{1}, y_{2}\right) \tilde{p}\left(y_{1} \mid x_{1}\right)=\tilde{p}\left(x_{1}\right) \tilde{p}\left(y_{2} \mid x_{1}\right) \tilde{p}\left(y_{1} \mid x_{1}\right) \end{array} $$

于是有：

$$
\tilde{p}\left(x_{1}, y_{1}\right) \tilde{p}\left(y_{2} \mid x_{1}\right)=\tilde{p}\left(x_{1}, y_{2}\right) \tilde{p}\left(y_{1} \mid x_{1}\right)
$$

则在 $ x=x_{1} $ 这条平行于 $ y $ 轴的直线上，如果使用条件分布 $ \tilde{p}\left(y \mid x_{1}\right) $ 作为直线上任意两点之间的转移概率, 则这两点之间的转移满足细致平稳条件。
同理: 考察 $ y $ 坐标相同的两个点 $ A\left(x_{1}, y_{1}\right), C\left(x_{2}, y_{1}\right) $, 有 $ \tilde{p}\left(x_{1}, y_{1}\right) \tilde{p}\left(x_{2} \mid y_{1}\right)=\tilde{p}\left(x_{2}, y_{1}\right) \tilde{p}\left(x_{1} \mid y_{1}\right) $ 。在$ y=y_{1} $ 这条平行于 $ x $ 轴的直线上，如果使用条件分布 $ \tilde{p}\left(x \mid y_{1}\right) $ 作为直线上任意两点之间的转移概率, 则这两点之间的转移满足细致平稳条件。
可以构造状态空间上任意两点之间的转移概率矩阵$ \mathbf{Q}: $ 对于任意两点 $ A=\left(x_{A}, y_{A}\right), B=\left(x_{B}, y_{B}\right), \quad $ 令从$ A $ 转移到 $ B $ 的概率为 $ Q(A \rightarrow B) $ :
- 如果 $ x_{A}=x_{B}=x^{*} $, 则 $ Q(A \rightarrow B)=\tilde{p}\left(y_{B} \mid x^{*}\right) $ 。
- 如果 $ y_{A}=y_{B}=y^{*} $, 则 $ Q(A \rightarrow B)=\tilde{p}\left(x_{B} \mid y^{*}\right) $ 。
- 否则 $ Q(A \rightarrow B)=0 $ 。
采用该转移矩阵$Q$，可以证明：对于状态空间中任意两点$A,B$ ，都满足细致平稳条件：

$$
\tilde{p}(A) Q(A \rightarrow B)=\tilde{p}(B) Q(B \rightarrow A)
$$

于是这个二维状态空间上的马尔可夫链将收敛到平稳分布 $ \tilde{p}(x, y) $, 这就是吉布斯采样的原理。

算法:

输入：目标分布 $ \tilde{p}(\overrightarrow{\mathbf{x}}) $, 其中 $ \overrightarrow{\mathbf{x}} \in \mathbb{R}^{n} $

输出: 采样出的一个状态序列 $ \left\{\overrightarrow{\mathbf{x}}_{0}, \overrightarrow{\mathbf{x}}_{1}, \cdots\right\} $

算法步骤：

初始化：随机初始化 $ \overrightarrow{\mathbf{x}}_{0}, t=0 $ 。
执行迭代，迭代步骤如下:
- 将 $ x_{1}, \cdots, x_{i-1}, x_{i+1}, \cdots, x_{n} $ 保持不变, 计算条件概率:
  $$
  \tilde{p}\left(x_{i} \mid x_{1}=x_{t, 1}, \cdots, x_{i-1}=x_{t, i-1}, x_{i+1}=x_{t, i+1}, \cdots, x_{n}=x_{t, n}\right)
  $$
  
  该条件概率就是状态转移概率 $ Q(A \rightarrow B) $
- 根据条件概率 $ \tilde{p}\left(x_{i} \mid x_{1}=x_{t, 1}, \cdots, x_{i-1}=x_{t, i-1}, x_{i+1}=x_{t, i+1}, \cdots, x_{n}=x_{t, n}\right) $ 对分量$ x_{i} $ 进行采样，假设采样结果为 $ x_{t, i}^{*} $，则获得新样本:
  $$ \overrightarrow{\mathbf{x}}_{t+1}=\left(x_{t, 1}, \cdots, x_{t, i-1}, x_{t, i}^{*}, x_{t, i+1}, \cdots, x_{t, n}\right)^{T} $$
- 令 $ t \leftarrow t+1 $, 继续遍历。
最终返回一个状态序列 $ \left\{\overrightarrow{\mathbf{x}}_{0}, \overrightarrow{\mathbf{x}}_{1}, \cdots\right\} $ 。