Markov Chain Monte Carlo 采样算法

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

简介

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

——百度百科

背景

中的概率可以计算为：

如果函数

如果 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})] ​ 的无偏估计:

如果概率密度函数 \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) , 使其满足:

• 可以取 \alpha(i, j)=\tilde{p}(j) Q_{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 ​, 则有:

• 于是有:

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

将 \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) ​​ :

• 根据均匀分布从 (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)

• 考虑二维的情形：假设有概率分布 \tilde{p}(x, y) ​, 考察状态空间上 x ​ 坐标相同的两个点 A\left(x_{1}, y_{1}\right), B\left(x_{1}, y_{2}\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}(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 。
• 保持不变, 计算条件概率:

该条件概率就是状态转移概率

• 最终返回一个状态序列 \left\{\overrightarrow{\mathbf{x}}_{0}, \overrightarrow{\mathbf{x}}_{1}, \cdots\right\} 。

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

参考资料

• http://www.huaxiaozhuan.com/数学基础/chapters/4_monte_carlo.html
• https://www.cnblogs.com/pinard/p/6625739.html
• https://baike.baidu.com/item/马尔科夫链蒙特卡洛方法/20869120?fr=aladdin