马尔可夫链蒙特卡洛（MCMC）算法

在之前的推送中我们了解到什么是马尔可夫链（Markov Chain)。下面我们来介绍一下马尔可夫链蒙特卡洛算法（Markov Chain Monte Carlo)， 在此之前，我们需要回顾一下马尔可夫链的极限分布（limiting behavior)。

对于一个不可约非周期性的马尔可夫链，其转移矩阵为P，当经过t->inf 步之后，其状态概率收敛于固定值， 即：

则转移矩阵

我们记向量\pi=（\pi_{0},\pi_{1},...） 此时pi满足\pi=\pi P且是唯一解。

以下我们所提到的两种算法都用到马尔可夫链的极限分布。

马尔可夫链蒙特卡洛（MCMC）算法的产生是为了解决计算机产生随机数的问题。产生的随机数要服从一定的概率分布P(X)，当这个目标概率分布不太复杂时，比如均匀分布，计算机可以根据算法产生较好的伪随机数。比较著名的有线性同余随机数生成器（Linear congruential generator, 更加高级的有 Mersenne Twister。但是当这个分布比较复杂，比如涉及到多元随机变量，此时可能有人会想到当随机变量之间相互独立，P(X1, X2...)可以写成几个概率分布相乘的形式。但是在实际生活中，随机变量之间一般是有联系的，此时我们就需要引入MCMC的两种算法： Metropolis-Hastings采样法和Gibbs采样法。

1 Metropolis-Hastings 采样

假定我们需要对目标分布p(x)进行取样。Metropolis-Hastings（M-H）算法的主要思路是构建一个马尔可夫链，其最终收敛的平稳分布恰好是我们想要的目标分布p(x)。

为了进一步讨论MH的收敛， 我们首先介绍一个概念：细致平稳条件（detailed balance condition)。细致平稳条件是指: 对于所有的 i, j

也就是说，在t达到平稳后，从t时间的状态向量pi转移到t+1时间的状态向量pi时， 对于任意状态i和j，（也就是pi向量的第i和j个元素）， 从i转移到j的量恰好等于从j转移到i的量。 这样就不难理解， 第i个状态的量就不会改变， 因为i是任意的，我们也可以说从t时间到t+1时间，pi向量的任意元素的量不发生改变。此时pi自然就达到平稳状态。这个条件只是充分不必要条件。

但是问题来了，如何能找到这样的一个转移矩阵P，能够最终满足细致平稳条件呢？通常情况下，P是不满足的，我们需要对其进行改造。假定对于任意一个转移矩阵， 我们定义它为Q，一般情况下：

算法：

给定现有状态Xt

1 采样Y~q(Xt,y)

2 采样U~U(0,1), 使得：

其中

3 重复以上步骤，当t足够大时，X_t近似服从目标分布

python代码：(2.7)

2 Gibbs 采样

当从条件分布采样比从联合分布采样更容易时， 我们常用Gibbs算法。Gibbs采样多用于多维度（2维以上），它的核心思想是：选定一维度i，固定其他维度的值不变对维度i进行采样，然后重复对所有维度做此处理。

其算法为：

1 给定X_t, 采集样本Y如下：

a. 从条件分布

中对Y1采样

b. 类似的，从条件分布

中对 Yi采样

c. 从条件分布

中对Yn采样。 本轮采样结束

2 令X_t+1=Y

Gibbs采样满足以上所说的细致平稳条件。以二维情形为例，假设坐标x恒定，y发生变化则

因此

同理，当我们固定y坐标，x发生变化，则

所以细致平稳条件成立。

当维度的选定是随机的时候，Gibbs采样可以看作是M-H算法的一个特例（接受率alpha＝1）。此时转移方程q为：

其中y=(x1,x2,x3...,xn), 而

因此

下面是一个简单的对二元正态分布采样的例子： （来自The Clever Machine）

对于一个二元正态分布：

其中

在Gibbs采样过程中，每一次迭代我们都先在x2的最新状态为条件下去对x1采样，然后反过来在刚采样得到的最新的x1为条件下对x2采样。其条件概率如下（数学部分省略）：

最终结果如下。

参考文献：

Simulation and the Monte Carlo Method (2nd edition) R.Y.Rubinstein, D.P.Kroese

The Clever Machine. https://theclevermachine.files.wordpress.com/2012/11/gibbssampler-2dnormal1.png