本文最后更新于: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 $, 则有:
- 于是有:
- 因此,即使提高了接受率, 细致平稳条件仍然成立。
将 $ \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) $ :
- 根据均匀分布从 $ (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) $,可以证明有:
- 于是有:
$$
\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 算法的特例, 它也使用马尔可夫链获取样本。
参考资料
- 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
文章链接:
https://www.zywvvd.com/notes/study/math/mcmc-alg/mcmc-alg/
“觉得不错的话,给点打赏吧 ୧(๑•̀⌄•́๑)૭”
微信支付
支付宝支付