Markov Chain Monte Carlo 采样算法

本文最后更新于:2021年11月1日 中午

作为一种随机采样方法,马尔科夫链蒙特卡罗(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 $, 此时满足细致平稳条件。

算法

输入:

  1. 先验转移概率矩阵$Q$

  2. 目标分布 $ \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\} $ 。

吉布斯采样 Gibbs sampling 有时被视作 MH 算法的特例, 它也使用马尔可夫链获取样本。

参考资料