MCMC算法

2022-05-01

字数统计: 7.8k | 阅读时长≈ 32 分钟

马尔科夫链蒙特卡罗方法（Markov Chain Monte Carlo，以下简称MCMC），是用于得到服从某一特定分布的随机采样的方法。很多复杂算法求解的基础，都使用到了 MCMC 算法。

本文内容主要参考张勤老师的《动物遗传育种中的计算方法》和刘建平的博客MCMC(一)蒙特卡罗方法

随机数产生方法

服从某一特定分布的随机变量的抽样值称为随机数，我们通常采用数学方法来产生随机数。

这种方法是利用一定的数学递推公式来产生随机数，即在给定一个任意的初值后，就可由此递推公式产生任意多的随机数，这种方法非常容易在计算机上实现。这种方法最大的缺陷在于给定初值后，以后所有的随机数便被唯一地确定下来了，而且这些随机数还存在周期现象，即随机序列达到一定长度后会出现重复，因此严格来说这些随机数并不是真正的随机数，故一般称之为伪随机数 (pseudo random number)。但由于初值可随机确定，因此其后的一系列数也可看作是随机的。

[0,1]均匀随机数的产生

[0,1]均匀随机数是在[0,1]区间上均匀分布随机变量的抽样值，其是产生其他一切分布随机数的基础。也就是说，任何其他分布的随机数都可通过对[0,1]均匀随机数进行某种转换得到。

目前产生[0,1]均匀随机数的最广泛的方法是线性同余法，也称为线性同余发生器 (linear congruential generator) 。

线性同余法

Lehmer (1951) 提出这种方法，它利用数论中的同余运算原理产生随机数，其递推公式为 \[ \begin{aligned} x_{i} &= (ax_{i-1} +c)(\mathrm{MOD} \quad m)\\ r_{i} &= x_{i}/m, \quad i=1,2,\cdots \\ \end{aligned} \]

其中 $a,c,m$ 是事先给定的参数，均为非负整数；$x_{i}$ 是在 [0,m] 区间取值的一个随机整数；$r_{i}$ 是在 [0,1] 区间取值的一个随机实数。

在给定一个任意非负整数的初值 $x_{0}$ 后，由上式求出整数序列 $x_{1},x_{2},\cdots$ 和实数序列 $r_{1},r_{2},\cdots$ 。

若 $c=0$ ，则称为乘同余法；若 $c>0$ ，则称为混同余法。

其他分布随机数的产生

任何其他分布的随机数都可以通过[0,1]均匀进行某种转换得到，例如正态分布可以通过著名的 Box-Muller 变换得到，根据下面公式通过两个[0,1]均匀随机数，我们可以得到两个彼此独立的标准正态分布随机数（缺证明，下式中的 $X_{1}$ 和 $X_{2}$ 为两个独立的[0,1]均匀随机数）。 \[ \begin{aligned} &Z_{1}=\sqrt{-2 \ln X_{1}} \cos \left(2 \pi X_{2}\right) \\ &Z_{2}=\sqrt{-2 \ln X_{1}} \sin \left(2 \pi X_{2}\right) \end{aligned} \] 下面介绍两种基本方法，反函数法与接受-拒绝采样法。

反函数法

假设我们已知某个分布的累积分布函数 $F(x)$ ，那么抽样值可以通过下式得到 \[ x = F^{-1}(u) \] 其中，$u$ 为[0,1]均匀随机数。

证明：因为 $u$ 服从[0,1]均匀分布，所以 $P(u \leq t) =t$ ，因此，我们有 \[ P\left(F^{-1}(u) \leqslant x\right)=\mathrm{P}\left(F\left(F^{-1}(u)\right) \leqslant F(x)\right)=\mathrm{P}(u \leqslant F(x))=\mathrm{F}(x) \] 所以 $F^{-1}(u)$ 的累积分布函数就是 $F(x)$ ，因此 $x = F^{-1}(u)$ 就是 $F(x)$ 的一个抽样值。

这种方法非常简单，但是要求随机变量累积分布函数的反函数存在，而且容易计算。在很多情况下，累积分布函数的反函数不存在，或者虽然存在但其计算却非常困难，此时这种方法并不适用。

接受-拒绝采样法

对于概率分布不常见的分布，假设我们已知其概率密度分布为 $p(x)$ ，此时我们用一个程序可采样的分布 $q(x)$ ，比如正态分布。具体过程如下：

首先我们找到一个常数 $k$ , 使得 $p(x) \leq kq(x)$ 恒成立，即 $p(x)$ 恒在 $kq(x)$ 下方。
采样得到 $q(x)$ 的一个样本 $x_{0}$ ，然后从均匀分布 $[0,kq(x_{0})]$ 中采样得到一个值 $u$ 。如果 $u > p(x_{0})$ ，则拒绝这一次抽样，否则则接受这次抽样。
重复以上过程，直到得到 $n$ 个抽样值。

整个过程从直观理解来看，就是从 $kq(x)$ 曲线下方的二维图形中随机抽取一个点，如果这个点正好落在 $p(x)$ 的曲线下方，则接受这个抽样，否则就拒绝这个抽样。

因此使用接受-拒绝采样法时，只有当 $kq(x)$ 与 $p(x)$ 曲线图形越接近时，抽样接受率越高，抽样效率越高。

蒙特卡罗方法

蒙特卡罗方法 (Monte Carlo) 是指使用随机数来解决很多计算问题的方法。

总的来说，蒙特卡罗方法处理的问题可以分为两大类：

第一类是确定性的数学问题。如计算定积分、解线性和非线性方程组，解积分方程和某些偏积分方程等。在这类问题中，因不含时间因素，故常将相应的模拟模型称为静态模型。

第二类是随机性问题，例如中子在介质中的扩散问题，运筹学中的库存问题等。

用蒙特卡罗方法计算定积分

用蒙特卡罗方法计算定积分是其最典型的应用，这里有两种算法，一是随机投点法，二是平均值法。下面以单重积分为例，分别进行介绍。

随机投点法

假设所求定积分为 \[ I= \int_{0}^{1} g(x) dx, \quad 0 \leq g(x) \leq 1 \] 因此，如果我们向 (0,0) (0,1) (1,1) (1,0) 这四个点构成的单位正方形均匀地投点，落在曲线 $g(x)$ 下方的概率就是我们要求的定积分，证明如下： \[ P(y \leq g(x)) = \int_{0}^{1} \int_{0}^{g(x)} dydx = \int_{0}^{1} g(x) dx \] 因此，我们可以按照以下步骤计算这个积分：

独立地产生两个 [0,1]均匀随机数 $x_{i}$ 和 $y_{i}$ ，则 $(x_{i}, y_{i})$ 就构成了在单位正方形中地一个均匀分布的随机点
检验该随机点是否落在曲线 $g(x)$ 下方

重复上面两个步骤 $N$ 次，设有 $m$ 次随机点落在曲线 $g(x)$ 下方，则当 $N$ 足够大时，我们有 \[ I \approx m/N \] 我们可以将上面的方法一般化，假设我们所求的定积分为 \[ I= \int_{a}^{b} g(x) dx, \quad L \leq g(x) \leq M \] 我们需要先对这个积分进行变换。由 $L \leq g(x) \leq M$ 可得 $0 \leqslant \frac{g(x)-L}{M-L} \leqslant 1$ ，因此我们有 \[ \begin{aligned} I &=\int_{a}^{b} g(x) \mathrm{d} x \\ &=\int_{a}^{b}(M-L) \frac{g(x)-L+L}{M-L} \mathrm{~d} x \\ &=(M-L) \int_{a}^{b} \frac{g(x)-L}{M-L} \mathrm{~d} x+L(b-a) \end{aligned} \] 令 $x = a+(b-a)z$ ，有 \[ \begin{aligned} &I=(M-L) \int_{0}^{1} \frac{g(a+(b-a) z)-L}{M-L}(b-a) d z+L(b-a) \\ & \\ &=(M-L)(b-a) \int_{0}^{1} g^{*}(z) d z+L(b-a) \end{aligned} \] 其中 $g^{*}(z) = \frac{g(a+(b-a) z)-L}{M-L}$ ，它满足 $0 \leq g^{*}(x) \leq 1$ ，因此积分 $ _{0}^{1} g^{*}(z) d z$ 可以按照上面的方法进行计算。

平均值法

假设我们所求的定积分为 \[ I= \int_{a}^{b} g(x) dx \] 设 $R$ 是 [a,b] 区间上的均匀分布随机变量，则 \[ E(g(R))=\int_{a}^{b} g(x) p(x) \mathrm{d} x=\int_{a}^{b} g(x) \frac{1}{b-a} \mathrm{~d} x=\frac{1}{b-a} I \] 于是，我们有 \[ I =(b-a)E(g(R)) \] 因此，我们可以用以下方法计算积分的近似解

产生 [a,b] 区间上的均匀分布随机数 $r_{i}$ ，$i=1,2,\cdots,N$
计算 $g(r_{i})$

当 $N$ 足够大时，我们有 \[ I \approx \frac{b-a}{N} \sum_{i=1}^{N} g(r_{i}) \] 其实这种做法就是在 [a,b] 区间上随机采样 $N$ 个点，用这 $N$ 个点的函数值均值代表 [a,b] 区间上的所有函数值，即用一个矩形的面积代表原来的曲线下方的面积。

这里我们有一个隐藏的假设是 $x$ 在 [a,b] 区间是均匀分布的，而大多数时候这个假设均不成立。此时我们进一步假设 $x$ 在 [a,b] 区间上的概率密度函数为 $p(x)$ ，此时我们有蒙特卡罗方法的一般形式： \[ I=\int_{a}^{b} g(x) d x=\int_{a}^{b} \frac{g(x)}{p(x)} p(x) d x = E(\frac{g(x)}{p(x)}) \approx \frac{1}{N} \sum_{i=1}^{N} \frac{g\left(x_{i}\right)}{p\left(x_{i}\right)} \] 因此，平均值法计算积分的一般过程如下：

从概率密度函数为 $p(x)$ 的分布中生成随机数 $ x_{i}$ ，$i=1,2,\cdots,N$
计算 $g\left(x_{i}\right)/p\left(x_{i}\right)$

当 $N$ 足够大时，我们有 \[ I \approx \frac{1}{N} \sum_{i=1}^{N} \frac{g\left(x_{i}\right)}{p\left(x_{i}\right)} \] 可以看出当我们假设$x$ 在 [a,b] 区间满足均匀分布时，我们将 $p(x) = \frac{1}{b-a}$ 带入上式，得到 \[ \begin{aligned} \frac{1}{N} \sum_{i=1}^{N} \frac{g\left(x_{i}\right)}{p\left(x_{i}\right)} &= \frac{1}{N} \sum_{i=1}^{N} \frac{g\left(x_{i}\right)}{1/(b-a)} \\ &= \frac{b-a}{N} \sum_{i=1}^{N} g(x_{i}) \end{aligned} \] 我们得到了均匀分布的估计公式，即均匀分布是上式的一个特列。

小结

蒙特卡罗方法的核心在于生成特定分布的随机数，常用的分布可以通过对于 [0,1] 均匀随机数进行转换得到。通过采用接受-拒绝采样，我们也可以得到一些不常见分布的随机数。

但是在很多时候，有一些分布我们还是很难得到其随机数，比如：

1）对于一些二维分布$p(x,y)$，这个函数本身计算困难，但是我们能得到其条件分布 $p(x|y)$ 和 $p(y|x)$ ，这时我们无法用接受-拒绝采样得到其样本集。

2）对于一些高维的分布 $p(x_{1},x_{2},...,x_{n})$ ，我们要找到一个合适的 $q(x)$ 和 $k$ 非常困难。

此时我们需要马尔科夫链的方法帮助我们进行采样。

马尔科夫链

马尔科夫链 (Markov chain) 是一种最简单，应用最多的随机过程。随机过程是依赖于时间的一组随机变量。马尔科夫链要求具备“无记忆”的性质：下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关，表示为 \[ P\left(X_{t+1}=x \mid X_{t}, X_{t-1}, \cdots\right)=P\left(X_{t+1}=x \mid X_{t}\right) \] 我们将系统随着时间的变化而发生的状态改变称为系统状态的转移，将已知在时刻 $t$ 系统处于状态 $s_{i}$ 的条件下，在下一时刻处于状态 $s_{j}$ 的概率称为转移概率，表示为 $P_{ij}(t, t+1)$ ，即 \[ P_{ij}(t, t+1) = P\{ X(t+1) = s_{j} | X(t) = s_{i}\} \] 如果系统由状态 $i$ 在下一时刻转移为状态 $j$ 的概率与系统当前所处的时刻无关，则称其为齐次的 (homogeneous) 。对于齐次的马尔科夫链，我们可以将转移概率简单地表示为 $P_{ij}$ 。一般来说，我们涉及到的马尔科夫链均为齐次的，因此下文均默认为齐次的马尔科夫链。

此时我们将所有的转移概率放到一个矩阵中，称为转移概率矩阵，如下。 \[ \boldsymbol{P}=\left[\begin{array}{cccc} P_{11} & P_{12} & P_{13} & \cdots \\ P_{21} & P_{22} & P_{23} & \cdots \\ P_{31} & P_{32} & P_{33} & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ \end{array}\right] \] 转移概率矩阵的每一行之和为1，即 \[ \sum_{j} P_{ij} = 1 \] 举个例子，假设天气变化是一个随机过程，假设天气的状态空间为 (有雨，晴天，多云)，再假设天气变化是一个马尔科夫链过程，即明天出现何种天气的概率仅与今天的天气有关。设状态转移矩阵为 \[ \boldsymbol{P}=\left[\begin{array}{ccc} 0.5 & 0.25 & 0.25 \\ 0.5 & 0 & 0.5 \\ 0.25 & 0.25 & 0.5 \end{array}\right] \] 其中，第一行表示如果今天有雨，则明天为有雨，晴天，多云的概率分别为 0.5，0.25 和 0.25 。

假设系统初始的概率分布向量为 $\boldsymbol{\pi}_{0}=\left[\pi_{0}(1), \pi_{0}(2), \pi_{0}(3)\right]$ ，那么第一天的天气概率分布向量为 $\boldsymbol{\pi}_{1} = \boldsymbol{\pi}_{0} \mathbf{P}$ ，第二天的分布向量为 $\boldsymbol{\pi}_{2} = \boldsymbol{\pi}_{0} \mathbf{P}^{2}$，依次类推，第 $n$ 天的分布向量为 $\boldsymbol{\pi}_{n} = \boldsymbol{\pi}_{0} \mathbf{P}^{n}$ .

我们发现随着天数的增加，$\boldsymbol{\pi}_{n}$ 会逐渐变成一个稳定不变的分布，此时我们称马尔科夫链收敛，收敛时系统状态空间的概率分布称为平稳分布 (stationary distribution)，表示为 $\boldsymbol{\pi}$ 。

在本例中为 $\boldsymbol{\pi} = (0.4 ,0.2, 0.4)$ ，而且与初始值无关。进一步发现其原因是由于 $\mathbf{P}^{n}$ 在 $n =7$ 之后就变成了一个稳定不变的矩阵，即 \[ \mathbf{p}^{7}=\left[\begin{array}{lll} 0.4 & 0.2 & 0.4 \\ 0.4 & 0.2 & 0.4 \\ 0.4 & 0.2 & 0.4 \end{array}\right]=\left[\begin{array}{l} \boldsymbol{\pi} \\ \boldsymbol{\pi} \\ \boldsymbol{\pi} \end{array}\right] \] 但是不是所有的马尔科夫链都能达到平稳分布。

性质

不可约：如果从马尔科夫链中的任一状态出发，都能以大于零的概率经过有限次转移变为其他任一状态，则称该马尔科夫链是不可约的 (irreducible) 。假设 $P_{ij}(k)$ 为从状态 $i$ 经过 k 次 ($k\geq 1$) 转移变成状态 $j$ 的概率，不可约是说对于任意的 $i$ 和 $j$ ，必定存在一个正整数 $k$ ，使得 $P_{ij}(k) > 0$ 。换句话说，如果马尔科夫链是不可约的，那么它的任何两个状态是连通的（两个状态彼此均可以通过有限次转移变为另一种状态，记为 $s_{i} \leftrightarrow s_{j}$ ），即可以经过一次或多次转移从某一种状态转移到任意一种状态。

非周期：如果所有满足 $P_{ii}(k) > 0$ 的 $k$ 中没有大于1的公因子，则称状态 $i$ 是非周期的 (aperiodic)。举个例子，假设存在一个马尔科夫链，其状态空间为 (1,2,3)，而状态之间的转移只能是 1-2-3-1-2-3 ，即系统状态每经过三次转移就会发生重复，假设我们从状态 2 开始，则经过 3,6,9… 次转移之后会重复出现状态 2 ，因此我们说状态 2 是周期性的，其周期为3。由周期性的定义存在一下推论：

推论：若状态A与状态B连通，则A与B周期相同。
推论：若不可约的马尔科夫链有周期性状态A，则该马尔科夫链的所有状态为周期性状态

只有不可约和非周期的马尔科夫链才能收敛于平衡状态，并且存在以下性质

\[ \lim _{n \rightarrow \infty} P^{n}=\left(\begin{array}{ccccc} \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \pi(1) & \pi(2) & \ldots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \pi(1) & \pi(2) & \ldots & \pi(j) & \cdots \\ \cdots & \ldots & \cdots & \cdots & \ldots \end{array}\right) \]
平稳分布 $\boldsymbol{\pi}$ 是 $\boldsymbol{\pi} \mathbf{P} = \boldsymbol{\pi}$ 的唯一非负解。

基于马尔科夫链采样

假如我们想得到某个马尔科夫链的平稳分布的随机数，如果我们得到了其转移概率矩阵，那么我们就可以轻松做到这一点。

我们可以从任意一个初始概率分布 $\boldsymbol{\pi}_{0}$ ，假设我们认为经过 $n$ 轮后马尔科夫链收敛于平稳分布，那么之后采样得到的随机数就是我们需要服从平稳分布的随机数，其过程如下：

输入转移概率矩阵 $\mathbf{P}$ ，设定状态转移次数阈值为 $n_{1}$ ，需要的样本数为 $n_{2}$
从任意简单概率分布采样得到初始状态值 $x_{0}$
从条件概率分布 $P(x|x_{t})$ 采用得到 $x_{t+1}$ ，$t = 0,1,\cdots, n_{1}+n_{2}-1$ 。

最终，样本集合 ($x_{n_{1}+1},x_{n_{1}+2},\cdots,x_{n_{1}+n_{2}}$ ) 即为我们需要的满足平稳分布的随机数。

但是随意给定一个平稳分布 $\boldsymbol{\pi}$ ，我们如何得到其相应的转移概率矩阵 $\mathbf{P}$ 呢？此时我们就需要将马尔科夫链和蒙特卡洛方法结合起来，即MCMC 方法。

MCMC 算法

在解决如果得到某个平稳分布相应的转移概率矩阵 $\mathbf{P}$ 这个问题之前，我们需要先看马尔科夫链的细致平稳条件

马尔科夫链的细致平稳条件

如果非周期马尔科夫链的转移概率矩阵 $\mathbf{P}$ 和概率分布 $\boldsymbol{\pi}(x)$ 对于所有的 $i,j$ ，均满足 \[ \pi(i) P(i, j)=\pi(j) P(j, i) \] 则 $\boldsymbol{\pi}(x)$ 就是转移概率矩阵 $\mathbf{P}$ 的平稳分布。

证明如下，由细致平稳条件有 \[ \sum_{i=1}^{\infty} \pi(i) P(i, j)=\sum_{i=1}^{\infty} \pi(j) P(j, i) \] 其中我们有 \[ \sum_{i=1}^{\infty} \pi(j) P(j, i)=\pi(j) \sum_{i=1}^{\infty} P(j, i)=\pi(j) \] 因此我们得到 \[ \sum_{i=1}^{\infty} \pi(i) P(i, j)= \pi(j) \] 即 \[ \boldsymbol{\pi} \mathbf{P} = \boldsymbol{\pi} \] 因此， $\boldsymbol{\pi}(x)$ 就是转移概率矩阵 $\mathbf{P}$ 的平稳分布，证明完毕。

MCMC 采样

一般情况下，我们的目标平稳分布 $\boldsymbol{\pi}(x)$ 和我们随机找的某一个马尔科夫链转移概率矩阵 $\mathbf{Q}$ 不满足细致平稳条件，即 \[ \pi(i) Q(i, j) \neq \pi(j) Q(j, i) \] 我们可以对上式做一个改造，使得细致平稳条件成立，方法式引入 $\alpha(i, j)$ 和 $\alpha(j, i)$ ，使得 \[ \pi(i) Q(i, j) \alpha(i, j)=\pi(j) Q(j, i) \alpha(j, i) \] 那么什么样的 $\alpha(i, j)$ 和 $\alpha(j, i)$ 可以使得上式成立呢？最简单的就是只要满足下面两个式子即可： \[ \begin{aligned} &\alpha(i, j)=\pi(j) Q(j, i) \\ &\alpha(j, i)=\pi(i) Q(i, j) \end{aligned} \] 这样，我们就得到了我们的目标平稳分布 $\boldsymbol{\pi}(x)$ 对应的转移概率矩阵 $\mathbf{P}$ ，其元素满足 \[ P(i, j)=Q(i, j) \alpha(i, j) \] 也就是说，我们的目标矩阵 $\mathbf{P}$ 可以通过任意一个转移概率矩阵 $\mathbf{Q}$ 构建得到，这里的 $\alpha(i, j)$ 我们一般称为接受率，这个思路很像上面的接受-拒绝采样。

MCMC 采样的步骤如下：

输入任意的转移概率矩阵 $\mathbf{Q}$ ，平稳分布 $\boldsymbol{\pi}(x)$ ，设定状态转移次数阈值 $n_{1}$ ，需要的样本个数 $n_{2}$
从任意简单概率分布采样得到初始状态值 $x_{0}$
$\text{for }t = 0,1,\cdots, n_{1}+n_{2}-1$ 。
1. 从条件概率分布 $Q(x|x_{t})$ 采用得到 $x_{*}$
2. 从[0,1]均匀分布采样 $u$
3. 如果$u<\alpha\left(x_{t}, x_{*}\right)=\pi\left(x_{*}\right) Q\left(x_{*}, x_{t}\right)$ ，则接受转移，即 $x_{t+1}=x_{*}$ ；否则不接受转移，即 $x_{t+1}=x_{t}$

最终，样本集合 ($x_{n_{1}+1},x_{n_{1}+2},\cdots,x_{n_{1}+n_{2}}$ ) 即为我们需要的满足平稳分布的随机数。

但是这个采样算法很难在实践中应用，因为这里的接受率 $\alpha(i, j)$ 可能非常的小，这样导致大部分的采样值都被拒绝转移，采样效率很低。

我们采用M-H 采样方法可以提高采样效率。

M-H 采样

M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hasting扩展，其解决了接受率过低的问题。

我们回到MCMC采样的细致平稳条件： \[ \pi(i) Q(i, j) \alpha(i, j)=\pi(j) Q(j, i) \alpha(j, i) \] 我们的问题是接受率 $\alpha(i, j)$ 太小了，比如 $\alpha(i, j)=0.1$ ， $\alpha(j, i)=0.2$ ，此时如果等式的两边同时扩大 5 倍，此时接受率提高到了 0.5，但是细致平稳条件依然是满足的。

因此，我们的接受率可以改为 \[ \alpha(i, j)=\min \left\{\frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}, 1\right\} \] 其实这种做法可以换种方式解释，我们回到一开始我们说 $\mathbf{Q}$ 不满足细致平稳条件，即 \[ \pi(i) Q(i, j) \neq \pi(j) Q(j, i) \] 这可以分为两种情况 \[ \pi(i) Q(i, j) > \pi(j) Q(j, i) \]

\[ \pi(i) Q(i, j) < \pi(j) Q(j, i) \]

对于第一种情况，我们可以在左侧乘以概率 $\alpha(i, j)= \frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)} < 1$ ，使得细致平稳条件成立 \[ \pi(i) Q(i, j)\alpha(i, j) = \pi(j) Q(j, i) \] 对于第二种情况，我们可以在右侧乘以概率 $\alpha(j, i)= \frac{\pi(i) Q(i, j)}{\pi(j) Q(j, i)} < 1$ ，使得细致平稳条件成立 \[ \pi(i) Q(i, j) = \pi(j) Q(j, i)\alpha(j, i) \] 第一种情况的接受率为 $\alpha(i, j)= \frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}$ ，第二种情况的接受率为 $\alpha(i, j)= 1$ ，因此我们得到 \[ \alpha(i, j)=\min \left\{\frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}, 1\right\} \] 通过这个微小的改造，M-H采样算法过程如下

输入任意的转移概率矩阵 $\mathbf{Q}$ ，平稳分布 $\boldsymbol{\pi}(x)$ ，设定状态转移次数阈值 $n_{1}$ ，需要的样本个数 $n_{2}$
从任意简单概率分布采样得到初始状态值 $x_{0}$
$\text{for }t = 0,1,\cdots, n_{1}+n_{2}-1$ 。
1. 从条件概率分布 $Q(x|x_{t})$ 采用得到 $x_{*}$
2. 从[0,1]均匀分布采样 $u$
3. 如果$u<\alpha\left(x_{t}, x_{*}\right) = \min \left\{\frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}, 1\right\}$ ，则接受转移，即 $x_{t+1}=x_{*}$ ；否则不接受转移，即 $x_{t+1}=x_{t}$

最终，样本集合 ($x_{n_{1}+1},x_{n_{1}+2},\cdots,x_{n_{1}+n_{2}}$ ) 即为我们需要的满足平稳分布的随机数。

很多时候，我们选择的转移概率矩阵 $\mathbf{Q}$ 都是对称矩阵，即 $Q(i, j)=Q(j, i)$ ，此时我们的接受率可以进一步简化为 \[ \alpha(i, j)=\min \left\{\frac{\pi(j)}{\pi(i)}, 1\right\} \]

小结

M-H算法有下面几个问题：

对于高维度数据，由于每一次抽样均需要计算 $\frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}$ ，计算效率低
由于 $\alpha(i, j) < 1$，仍然有可能拒绝一部分抽样。
有时，对于高纬度数据，我们可能不知道其联合分布，但是可以方便计算各个特征之间的条件概率分布，此时无法使用 M-H 抽样。

吉布斯抽样 (Gibbs) 解决了上面的几个问题，其接受概率为1，即不拒绝任何采样，而且其抽样过程中只计算条件分布。

Gibbs 抽样

二维抽样

我们回到细致平稳条件，我们有 \[ \pi(i) P(i, j)=\pi(j) P(j, i) \] 我们假设我们的数据是二维数据，平稳分布为 $\pi\left(x_{1}, x_{2}\right)$ ，观测第一个特征相同的两个点 $A\left(x_{1}^{(1)}, x_{2}^{(1)}\right)$ 和 $B\left(x_{1}^{(1)}, x_{2}^{(2)}\right)$ ，根据概率论的知识，容易发现下面两个式子成立 \[ \begin{aligned} &\pi\left(x_{1}^{(1)}, x_{2}^{(1)}\right) \pi\left(x_{2}^{(2)} \mid x_{1}^{(1)}\right)=\pi\left(x_{1}^{(1)}\right) \pi\left(x_{2}^{(1)} \mid x_{1}^{(1)}\right) \pi\left(x_{2}^{(2)} \mid x_{1}^{(1)}\right) \\ &\pi\left(x_{1}^{(1)}, x_{2}^{(2)}\right) \pi\left(x_{2}^{(1)} \mid x_{1}^{(1)}\right)=\pi\left(x_{1}^{(1)}\right) \pi\left(x_{2}^{(2)} \mid x_{1}^{(1)}\right) \pi\left(x_{2}^{(1)} \mid x_{1}^{(1)}\right) \end{aligned} \] 由于两个等式的右手项相等，因此我们有 \[ \pi\left(x_{1}^{(1)}, x_{2}^{(1)}\right) \pi\left(x_{2}^{(2)} \mid x_{1}^{(1)}\right)=\pi\left(x_{1}^{(1)}, x_{2}^{(2)}\right) \pi\left(x_{2}^{(1)} \mid x_{1}^{(1)}\right) \] 也就是说 \[ \pi(A) \pi\left(x_{2}^{(2)} \mid x_{1}^{(1)}\right)=\pi(B) \pi\left(x_{2}^{(1)} \mid x_{1}^{(1)}\right) \] 联合上面的细致平稳条件，我们发现在 $ x_{1}=x_{1}^{(1)}$ 这条直线上，我们只需要选择 $P(A,B) = \pi\left(x_{2}^{(B)} \mid x_{1}^{(1)}\right)$ ，则任意两个点均满足细致平稳条件。

同样的道理，在在 $ x_{2}=x_{2}^{(1)}$ 这条直线上，我们只需要选择 $P(A,C) = \pi\left(x_{1}^{(C)} \mid x_{2}^{(1)}\right)$ ，则任意两个点也满足细致平稳条件。

基于上面的发现，我们可以构建 $\pi\left(x_{1}, x_{2}\right)$ 的转移概率矩阵为 \[ \begin{aligned} P(A \rightarrow B)=& \pi\left(x_{2}^{(B)} \mid x_{1}^{(1)}\right) \text { if } x_{1}^{(A)}=x_{1}^{(B)}=x_{1}^{(1)} \\ P(A \rightarrow C)=& \pi\left(x_{1}^{(C)} \mid x_{2}^{(1)}\right) \text { if } x_{2}^{(A)}=x_{2}^{(C)}=x_{2}^{(1)} \\ & P(A \rightarrow D)=0 \text { else } \end{aligned} \] 我们可以轻松验证对于二维平面上的任意两点 $E,F$ ，细致平稳条件均成立，即： \[ \pi(E) P(E \rightarrow F)=\pi(F) P(F \rightarrow E) \] 于是这个二维空间的马尔科夫链将收敛到平稳分布 $\pi\left(x_{1}, x_{2}\right)$

二维 Gibbs 抽样具体过程如下：

输入平稳分布 $\pi\left(x_{1}, x_{2}\right)$ ，设定状态转移次数阈值 $n_{1}$ ，需要的样本个数 $n_{2}$
从任意简单概率分布采样得到初始状态值 $x_{1}^{(0)}$ 和 $x_{2}^{(0)}$
$\text{for }t = 0,1,\cdots, n_{1}+n_{2}-1$ 。
1. 从条件概率分布 $P(x_{1}|x_{2}^{(t)})=\pi(x_{1}|x_{2}^{(t)})$ 采样得到 $x^{(t+1)}_{1}$
2. 从条件概率分布 $P(x_{2}|x_{1}^{(t+1)}) = \pi(x_{2}|x_{1}^{(t+1)})$ 采样得到 $x^{(t+1)}_{2}$

最终，样本集合 $\left\{\left(x_{1}^{\left(n_{1}+1\right)}, x_{2}^{\left(n_{1}+1\right)}\right),\left(x_{1}^{\left(n_{1}+2\right)}, x_{2}^{\left(n_{1}+2\right)}\right), \ldots,\left(x_{1}^{\left(n_{1}+n_{2}\right)}, x_{2}^{\left(n_{1}+n_{2}\right)}\right)\right\}$ 即为我们需要的满足平稳分布的随机数。

整个采样过程中，我们不停轮换坐标轴，即 \[ \left(x_{1}^{(1)}, x_{2}^{(1)}\right) \rightarrow\left(x_{1}^{(2)}, x_{2}^{(1)}\right) \rightarrow\left(x_{1}^{(2)}, x_{2}^{(2)}\right) \rightarrow \ldots \rightarrow\left(x_{1}^{\left(n_{1}+n_{2}\right)}, x_{2}^{\left(n_{1}+n_{2}\right)}\right) \]

多维抽样

我们从二维扩展到多维，假设我们有一个 $n$ 维的平稳分布 $\pi\left(x_{1}, x_{2}, \ldots x_{n}\right)$ ，我们可以通过对 $n$ 个坐标轴轮换采样，来得到新的样本。具体的抽样过程如下：

输入平稳分布 $\pi\left(x_{1}, x_{2}, \ldots, x_{n}\right)$ 或者对应的所有特征的条件概率分布，设定状态转移次数阈值 $n_{1}$ ，需要的样本个数 $n_{2}$
随机初始化初始状态值 $\left(x_{1}^{(0)}, x_{2}^{(0)}, \ldots, x_{n}^{(0)}\right)$
for $t=0$ to $n_{1}+n_{2}-1$ :
1. 从条件概率分布 $P\left(x_{1} \mid x_{2}^{(t)}, x_{3}^{(t)}, \ldots, x_{n}^{(t)}\right)$ 中采样得到样本 $x_{1}^{t+1}$
2. 从条件概率分布 $P\left(x_{2} \mid x_{1}^{(t+1)}, x_{3}^{(t)}, x_{4}^{(t)}, \ldots, x_{n}^{(t)}\right)$ 中采样得到样本 $x_{2}^{t+1}$
3. …
4. 从条件概率分布 $P\left(x_{j} \mid x_{1}^{(t+1)}, x_{2}^{(t+1)}, \ldots, x_{j-1}^{(t+1)}, x_{j+1}^{(t)} \ldots, x_{n}^{(t)}\right)$ 中采样得到样本 $x_{j}^{t+1}$
5. …
6. 从条件概率分布 $P\left(x_{n} \mid x_{1}^{(t+1)}, x_{2}^{(t+1)}, \ldots, x_{n-1}^{(t+1)}\right)$ 中采样得到样本 $x_{n}^{t+1}$

样本集 $\left\{\left(x_{1}^{\left(n_{1}+1\right)}, x_{2}^{\left(n_{1}+1\right)}, \ldots, x_{n}^{\left(n_{1}+1\right)}\right), \ldots,\left(x_{1}^{\left(n_{1}+n_{2}\right)}, x_{2}^{\left(n_{1}+n_{2}\right)}, \ldots, x_{n}^{\left(n_{1}+n_{2}\right)}\right)\right\}$ 即为我们需要的平稳分布对应的样本集。

小结

由于 Gibbs 抽样对于高维度特征具有优势，因此一般所说的 MCMC 采样都是用的 Gibbs 抽样。但是 Gibbs 抽样要求至少有两个维度，因此一维的数据采样只能采用 M-H 采样。

收敛性判断

在MCMC 算法中，无论是 M-H 算法还是 Gibbs 算法，我们均认为随着马尔科夫链的延长，它就越接近其平稳分布，我们称其为收敛。而收敛前的所有迭代是要被舍弃的，这些被舍弃的抽样称为 burn-in 期。但是我们如何判断一条链是否收敛呢？也就是说，上面的 $n_{1}$ 应该设定为多少了？

这似乎没有一个统一的答案，目前常用的有两种方法。

图式法

我们将 MCMC 链上的轨迹点进行可视化，横坐标是迭代次数，纵坐标是抽样值。如果经过一段时间的剧烈波动之后逐渐趋于平稳，则意味着这条链可能收敛了。

更为合理的做法是，用不同的参数初值构建两个甚至多个独立的链，观察不同链的轨迹点，随着链的延长，不同链开始发生重叠，甚至几乎完全重合，此时可以认为收敛。

需要说明的是，图式法不能保证所观察的链一定收敛了，为稳妥起见， burn-in 期稍微设大一点较好。

Geweke 检验

首先假定一个burn-in 期，在其后的抽样值按一定间隔取2个样本，如前 10% 和后 10% 的抽样值的平均数之间不存在显著差异，则可认为链收敛了。否则，需要延长burn-in 期，再重新进行检验。

参考文献

张勤. 动物遗传育种中的计算方法[M]. 科学出版社, 2007.
MCMC(一)蒙特卡罗方法
Brémaud P. Markov chains: Gibbs fields, Monte Carlo simulation, and queues[M]. Springer Science & Business Media, 2013.

版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！