马尔科夫链蒙特卡罗方法（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) 提出这种方法，它利用数论中的同余运算原理产生随机数，其递推公式为

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

在给定一个任意非负整数的初值 后，由上式求出整数序列 和实数序列 。

若 ，则称为乘同余法；若 ，则称为混同余法。

其他分布随机数的产生

任何其他分布的随机数都可以通过[0,1]均匀进行某种转换得到，例如正态分布可以通过著名的 Box-Muller 变换得到，根据下面公式通过两个[0,1]均匀随机数，我们可以得到两个彼此独立的标准正态分布随机数（缺证明，下式中的 和 为两个独立的[0,1]均匀随机数）。

下面介绍两种基本方法，反函数法与接受-拒绝采样法。

反函数法

假设我们已知某个分布的累积分布函数 ，那么抽样值可以通过下式得到

其中， 为[0,1]均匀随机数。

证明：因为 服从[0,1]均匀分布，所以 ，因此，我们有

所以 的累积分布函数就是 ，因此 就是 的一个抽样值。

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

接受-拒绝采样法

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

1. 首先我们找到一个常数 , 使得 恒成立，即 恒在 下方。

2. 采样得到 的一个样本 ，然后从均匀分布 中采样得到一个值 。如果 ，则拒绝这一次抽样，否则则接受这次抽样。

3. 重复以上过程，直到得到 个抽样值。

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

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

蒙特卡罗方法

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

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

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

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

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

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

随机投点法

假设所求定积分为

因此，如果我们向 (0,0) (0,1) (1,1) (1,0) 这四个点构成的单位正方形均匀地投点，落在曲线 下方的概率就是我们要求的定积分，证明如下：

因此，我们可以按照以下步骤计算这个积分：

1. 独立地产生两个 [0,1]均匀随机数 和 ，则 就构成了在单位正方形中地一个均匀分布的随机点
2. 检验该随机点是否落在曲线 下方

重复上面两个步骤 次，设有 次随机点落在曲线 下方，则当 足够大时，我们有

我们可以将上面的方法一般化，假设我们所求的定积分为

我们需要先对这个积分进行变换。由 可得 ，因此我们有

令 ，有

其中 ，它满足 ，因此积分 $ \int_{0}^{1} g^{*}(z) d z$ 可以按照上面的方法进行计算。

平均值法

假设我们所求的定积分为

设 是 [a,b] 区间上的均匀分布随机变量，则

于是，我们有

因此，我们可以用以下方法计算积分的近似解

1. 产生 [a,b] 区间上的均匀分布随机数 ，
2. 计算

当 足够大时，我们有

其实这种做法就是在 [a,b] 区间上随机采样 个点，用这 个点的函数值均值代表 [a,b] 区间上的所有函数值，即用一个矩形的面积代表原来的曲线下方的面积。

这里我们有一个隐藏的假设是 在 [a,b] 区间是均匀分布的，而大多数时候这个假设均不成立。此时我们进一步假设 在 [a,b] 区间上的概率密度函数为 ，此时我们有蒙特卡罗方法的一般形式 ：

因此，平均值法计算积分的一般过程如下：

1. 从概率密度函数为 的分布中生成随机数 $ x_{i}$ ，
2. 计算

当 足够大时，我们有

可以看出当我们假设 在 [a,b] 区间满足均匀分布时，我们将 带入上式，得到

我们得到了均匀分布的估计公式，即均匀分布是上式的一个特列。

小结

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

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

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

2）对于一些高维的分布 ，我们要找到一个合适的 和 非常困难。

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

马尔科夫链

马尔科夫链 (Markov chain) 是一种最简单，应用最多的随机过程。随机过程是依赖于时间的一组随机变量。马尔科夫链要求具备“无记忆”的性质：下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关，表示为

我们将系统随着时间的变化而发生的状态改变称为系统状态的转移，将已知在时刻 系统处于状态 的条件下，在下一时刻处于状态 的概率称为转移概率，表示为 ，即

如果系统由状态 在下一时刻转移为状态 的概率与系统当前所处的时刻无关，则称其为齐次的 (homogeneous) 。对于齐次的马尔科夫链，我们可以将转移概率简单地表示为 。一般来说，我们涉及到的马尔科夫链均为齐次的，因此下文均默认为齐次的马尔科夫链。

此时我们将所有的转移概率放到一个矩阵中，称为转移概率矩阵，如下。

转移概率矩阵的每一行之和为1，即

举个例子，假设天气变化是一个随机过程，假设天气的状态空间为 (有雨，晴天，多云)，再假设天气变化是一个马尔科夫链过程，即明天出现何种天气的概率仅与今天的天气有关。设状态转移矩阵为

其中，第一行表示如果今天有雨，则明天为有雨，晴天，多云的概率分别为 0.5，0.25 和 0.25 。

假设系统初始的概率分布向量为 ，那么第一天的天气概率分布向量为 ，第二天的分布向量为 ，依次类推，第 天的分布向量为 .

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

在本例中为 ，而且与初始值无关。进一步发现其原因是由于 在 之后就变成了一个稳定不变的矩阵，即

但是不是所有的马尔科夫链都能达到平稳分布。

性质

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

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

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

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

1. 

2. 平稳分布 是 的唯一非负解。

基于马尔科夫链采样

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

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

1. 输入转移概率矩阵 ，设定状态转移次数阈值为 ，需要的样本数为
2. 从任意简单概率分布采样得到初始状态值
3. 从条件概率分布 采用得到 ， 。

最终，样本集合 ( ) 即为我们需要的满足平稳分布的随机数。

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

MCMC 算法

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

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

如果非周期马尔科夫链的转移概率矩阵 和概率分布 对于所有的 ，均满足

则 就是转移概率矩阵 的平稳分布。

证明如下，由细致平稳条件有

其中我们有

因此我们得到

即

因此， 就是转移概率矩阵 的平稳分布，证明完毕。

MCMC 采样

一般情况下，我们的目标平稳分布 和我们随机找的某一个马尔科夫链转移概率矩阵 不满足细致平稳条件，即

我们可以对上式做一个改造，使得细致平稳条件成立，方法式引入 和 ，使得

那么什么样的 和 可以使得上式成立呢？最简单的就是只要满足下面两个式子即可：

这样，我们就得到了我们的目标平稳分布 对应的转移概率矩阵 ，其元素满足

也就是说，我们的目标矩阵 可以通过任意一个转移概率矩阵 构建得到，这里的 我们一般称为接受率，这个思路很像上面的接受-拒绝采样。

MCMC 采样的步骤如下：

1. 输入任意的转移概率矩阵 ，平稳分布 ，设定状态转移次数阈值 ，需要的样本个数
2. 从任意简单概率分布采样得到初始状态值
3. 。
   1. 从条件概率分布 采用得到
   2. 从[0,1]均匀分布采样
   3. 如果 ，则接受转移，即 ；否则不接受转移，即

最终，样本集合 ( ) 即为我们需要的满足平稳分布的随机数。

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

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

M-H 采样

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

我们回到MCMC采样的细致平稳条件：

我们的问题是接受率 太小了，比如 ， ，此时如果等式的两边同时扩大 5 倍，此时接受率提高到了 0.5，但是细致平稳条件依然是满足的。

因此，我们的接受率可以改为

其实这种做法可以换种方式解释，我们回到一开始我们说 不满足细致平稳条件，即

这可以分为两种情况

对于第一种情况，我们可以在左侧乘以概率 ，使得细致平稳条件成立

对于第二种情况，我们可以在右侧乘以概率 ，使得细致平稳条件成立

第一种情况的接受率为 ，第二种情况的接受率为 ，因此我们得到

通过这个微小的改造，M-H采样算法过程如下

1. 输入任意的转移概率矩阵 ，平稳分布 ，设定状态转移次数阈值 ，需要的样本个数
2. 从任意简单概率分布采样得到初始状态值
3. 。
   1. 从条件概率分布 采用得到
   2. 从[0,1]均匀分布采样
   3. 如果 ，则接受转移，即 ；否则不接受转移，即

最终，样本集合 ( ) 即为我们需要的满足平稳分布的随机数。

很多时候，我们选择的转移概率矩阵 都是对称矩阵，即 ，此时我们的接受率可以进一步简化为

小结

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

1. 对于高维度数据，由于每一次抽样均需要计算 ，计算效率低
2. 由于 ，仍然有可能拒绝一部分抽样。
3. 有时，对于高纬度数据，我们可能不知道其联合分布，但是可以方便计算各个特征之间的条件概率分布，此时无法使用 M-H 抽样。

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

Gibbs 抽样

二维抽样

我们回到细致平稳条件，我们有

我们假设我们的数据是二维数据，平稳分布为 ，观测第一个特征相同的两个点 和 ，根据概率论的知识，容易发现下面两个式子成立

由于两个等式的右手项相等，因此我们有

也就是说

联合上面的细致平稳条件，我们发现在 $ x_{1}=x_{1}^{(1)}$ 这条直线上，我们只需要选择 ，则任意两个点均满足细致平稳条件。

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

基于上面的发现，我们可以构建 的转移概率矩阵为

我们可以轻松验证对于二维平面上的任意两点 ，细致平稳条件均成立，即：

于是这个二维空间的马尔科夫链将收敛到平稳分布

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

1. 输入平稳分布 ，设定状态转移次数阈值 ，需要的样本个数
2. 从任意简单概率分布采样得到初始状态值 和
3. 。
   1. 从条件概率分布 采样得到
   2. 从条件概率分布 采样得到

最终，样本集合 即为我们需要的满足平稳分布的随机数。

整个采样过程中，我们不停轮换坐标轴，即

多维抽样

我们从二维扩展到多维，假设我们有一个 维的平稳分布 ，我们可以通过对 个坐标轴轮换采样，来得到新的样本。具体的抽样过程如下：

1. 输入平稳分布 或者对应的所有特征的条件概率分布，设定状态转移次数阈值 ，需要的样本个数
2. 随机初始化初始状态值
3. for to :
   1. 从条件概率分布 中采样得到样本
   2. 从条件概率分布 中采样得到样本
   3. …
   4. 从条件概率分布 中采样得到样本
   5. …
   6. 从条件概率分布 中采样得到样本

样本集 即为我们需要的平稳分布对应的样本集。

小结

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

收敛性判断

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

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

图式法

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

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

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

Geweke 检验

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

参考文献

1. 张勤. 动物遗传育种中的计算方法[M]. 科学出版社, 2007.

2. MCMC(一)蒙特卡罗方法

3. Brémaud P. Markov chains: Gibbs fields, Monte Carlo simulation, and queues[M]. Springer Science & Business Media, 2013.