3.5　马尔科夫链蒙特卡洛

在以贝叶斯方法为基础的各种机器学习技术中，常常需要计算后验概率，再通过最大后验概率（MAP）的方法进行参数推断和决策。然而，在很多时候，后验分布的形式可能非常复杂，这个时候寻找其中的最大后验估计或者对后验概率进行积分等计算往往非常困难，此时可以通过采样的方法来求解。这其中非常重要的一个基础就是马尔科夫链蒙特卡洛（Markov chain Monte Carlo，MCMC）。

3.5.1　重要性采样

前面已经详细介绍了基于随机算法进行定积分求解的技术。这里主要用到其中的平均值法。这里仅做简单回顾。

在计算定积分时，会把g（x）拆解成两项的乘积，即g（x）=f（x）p（x），其f（x）是某种形式的函数，而p（x）是关于随机变量X的概率分布（也就是PDF或PMF）。如此一来，上述定积分就可以表示为求f（x）的期望，即

当然，g（x）的分布函数可能具有很复杂的形式，仍然无法直接求解，这时就可以用采样的方法去近似。这时积分可以表示为

在贝叶斯分析中，蒙特卡洛积分运算常常被用来在对后验概率进行积分时做近似估计。比如要计算

便可以使用下面这个近似形式

不难发现，在利用蒙特卡洛法进行积分求解时，非常重要的一个环节就是从特定的分布中进行采样。这里的“采样”意思也就是生成满足某种分布的观测值。之前已经介绍过“逆采样”和“拒绝采样”等方法。

采样的目的很多时候都是为了做近似积分运算，前面的采样方法（逆采样和拒绝采样）都是先对分布进行采样，然后再用采样的结果近似计算积分。下面要介绍的另外一种方法“重要性采样”（Importance sampling）则两步并做一步，直接做近似计算积分。

现在的目标是计算下面的积分

E［f（x）］=∫f（x）p（x）dx

按照蒙特卡洛求定积分的方法，将从满足p（x）的概率分布中独立地采样出一系列随机变量x_i，然后便有

但是现在的困难是对满足p（x）的概率分布进行采样非常困难，毕竟实际中很多p（x）的形式都相当复杂。这时该怎么做呢？于是想到做等量变换，于是有

如果把其中的看成是一个新的函数h（x），则有

其中被称为是x_i的权重或重要性权重（Importance Weights）。所以这种采样的方法就被称为是重要性采样（Importance Sampling）。

如图3-19所示，在使用重要性采样时，并不会拒绝掉某些采样点，这与在使用拒绝采样时不同。此时，所有的采样点都将为我们所用，但是它们的权重是不同的。因为权重为，所以在图3-19中，当p（x_i）＞q（x_i）时，采样点x_i的权重就大（深色线在浅色线上方时），反之就小。重要性采样就是通过这种方式来从一个“参考分布”q（x）中获得“目标分布”p（x）的。

3.5.2　马尔科夫链蒙特卡洛的基本概念

马尔科夫链蒙特卡洛方法是一类用于从一个概率分布中进行采样的算法，该类方法以构造马尔科夫链为基础，而被构造的马尔科夫链分布应该同需要之分布等价。

MCMC构造马尔科夫链，使其稳态分布等于要采样的分布，这样就可以通过马尔科夫链来采样。这种等价如何理解是深入探讨具体操作方法之前，需要先攻克的一个问题。在此之前，希望读者对马尔科夫链已经有了一个比较清晰的认识。

现在，用下面的式子来表示每一步（时刻推进）中从状态s_i到状态s_j的转移概率

p（i，j）=p（i→j）=P（X_t+1=s_j|X_t=s_i）

这里的一步是指时间从时刻t过渡到下一时刻t+1。

马尔科夫链在时刻t处于状态j的概率可以表示为π_j（t）=P（X_t=s_j）。这里用向量π（t）来表示在时刻t各个状态的概率向量。在初始时刻，需要选取一个特定的π（0），通常情况下可以使向量中一个元素为1，其他元素均为0，即从某个特定的状态开始。随着时间的推移，状态会通过转移矩阵，向各个状态扩散。

马尔科夫链在t+1时刻处于状态s_i的概率可以通过t时刻的状态概率和转移概率来求得，并可通过查普曼-柯尔莫哥洛夫等式得到

用转移矩阵写成矩阵的形式如下

π（t+1）=π（t）P

其中，转移矩中的元素p（i，j）表示p（i→j）。因此，π（t）=π（0）P^t。此外，表示矩阵Pⁿ中第ij个元素，即。

一条马尔科夫链有一个平稳分布π^∗，是指给定任意一个初始分布π（0），经过有限步之后，最终都会收敛到平稳分布π^∗。平稳分布具有性质π^∗=π^∗P。可以结合本章前面谈到的矩阵极限的概念来理解马尔科夫链的平稳分布。

一条马尔科夫链拥有平稳分布的一个充分条件是对于任意两个状态i和j，其满足细致平衡（Detailed Balance）：p（j→i）π_j=p（i→j）π_i。可以看到，此时

所以π=πP，明显满足平稳分布的条件。如果一条马尔科夫链满足细致平衡，就说它是可逆的。

最后来总结一下MCMC的基本思想。在拒绝采样和重要性采样中，当前生成的样本点与之前生成的样本点之间是没有关系的，它的采样都是独立进行的。然而，MCMC是基于马尔科夫链进行的采样。这就表明，当前的样本点生成与上一时刻的样本点是有关的。如图3-20所示，假设当前时刻生成的样本点是x，下一次时刻采样到x′的（转移）概率就是p（x|x′），或者写成p（x→x′）。我们希望的是这个过程如此继续下去，最终的概率分布收敛到π（x），也就是要采样的分布。

易见，转移概率p（x|x′）与采样分布π（x）之间必然存在某种联系，因为希望依照p（x|x′）来进行转移采样最终能收敛到π（x）。那这个关系到底是怎么样的呢？首先，给定所有可能的x，然后求从x转移到x′的概率，所得之结果就是x′出现的概率，如果用公式表示即

π（x′）=∫π（x）p（x′|x）dx

这其实也就是查普曼-柯尔莫哥洛夫等式所揭示的关系。如果马尔科夫链是平稳的，那么基于该马尔科夫链的采样过程最终将收敛到平稳状态

π^∗（x′）=∫π^∗（x′）p（x|x′）dx′

而平稳状态的充分条件又是满足细致平衡

π（x）p（x′|x）dx=π（x）p（x|x′）

可见细致平衡的意思是说从x转移到x′的概率应该等于从x′转移到x的概率，在这样的情况下，采样过程将最终收敛到目标分布。也就是说，设计的MCMC算法，只要验证其满足细致平衡的条件，那么用这样的方法做基于马尔科夫链的采样，最终就能收敛到一个稳定状态，即目标分布。

实际应用中，有两个马尔科夫链蒙特卡洛采样算法十分常用，它们是Metropolis-Hastings（梅特罗波利斯-黑斯廷斯）算法与Gibbs Sampling采样（吉布斯），可以证明吉布斯采样是梅特罗波利斯-黑斯廷斯算法的一种特殊情况，而且两者都满足细致平衡的条件。

3.5.3　Metropolis-Hastings算法

对给定的概率分布，如果从中直接采样比较困难，那么Metropolis-Hastings算法就是一个从中获取一系列随机样本的MCMC方法。这里所说的Metropolis-Hastings算法和模拟退火方法中所使用的Metropolis-Hastings算法本质上是一回事。用于模拟退火和MCMC的原始算法最初是由Metropolis提出的，后来Hastings对其进行了推广。Hastings提出没有必要使用一个对称的参考分布（proposal distribution）。当然达到均衡分布的速度与参考分布的选择有关。

Metropolis-Hastings算法的执行步骤如下。

Metropolis-Hastings算法的执行步骤是先随便指定一个初始的样本点x⁰，u是来自均匀分布的一个随机数，x^∗是来自p分布的样本点，样本点是否从前一个位置跳转到x^∗，由α（x^∗）和u的大小关系决定。如果接受，则跳转，即新的采样点为x^∗，否则就拒绝，即新的采用点仍为上一轮的采样点。这个算法看起来非常简单，难免让人疑问照此步骤来执行，最终采样分布能收敛到目标分布吗？根据之前的讲解，可知要想验证这一点，就需要检查细致平衡是否满足。下面是具体的证明过程，不再赘言。

既然已经从理论上证明Metropolis-Hastings算法的确实可行，下面举一个简单的例子来看看实际中Metropolis-Hastings算法的效果。稍微有点不一样的地方是，这里并未出现π（x）p（x′|x）这样的后验概率形式，作为一个简单的示例，下面只演示从柯西分布中进行采样的做法。

柯西分布也叫作柯西-洛伦兹分布，它是以奥古斯丁·路易·柯西与亨德里克·洛伦兹名字命名的连续概率分布，其概率密度函数为

其中，x₀是定义分布峰值位置的位置参数，γ是最大值一半处的一半宽度的尺度参数。x₀=0且γ=1的特例称为标准柯西分布，其概率密度函数为

下面是在R中进行基于Metropolis-Hastings算法对柯西分布进行采样的示例代码。

图3-21所示为每次采样点的数值分布情况。

为了更清晰明确地看到采样结果符合预期分布，这里也绘制了分布的密度图并与实际的分布密度进行对照，如图3-22所示，深色实线是采样分布的密度图，而浅色虚线则是实际柯西分布的密度图，可见两者吻合地相当好。

当然，作为一个简单例子的开始，上面的例子中并没有涉及p（y|x），接下来的这个例子则会演示涉及后验概率的基于Metropolis-Hastings算法的MCMC。这里将用Metropolis-Hastings算法从瑞利分布（Rayleigh Distribution）中采样，瑞利分布的密度为

然后取自由度为X_t的卡方分布为参考分布。实现代码如下。

然后要验证一下，生产的采样数据是否真的符合瑞利分布。注意在R中使用瑞利分布的相关函数，需要加装VGAM包。下面的代码可以让读者直观地感受到采样结果的分布情况。

从图3-23中不难看出，采样点分布确实符合预期。当然，这仅仅是一个演示用的小例子。显然，它并不高效。因为采用自由度为X_t的卡方分布来作为参考函数，大约有40%的采样点都被拒绝掉了。如果换做其他更加合适的参考函数，可以大大提升采样的效率。

3.5.4　Gibbs采样

Gibbs采样是一种用以获取一系列观察值的MCMC算法，这些观察值近似于来自一个指定的多变量概率分布，但从该分布中直接采样较为困难。Gibbs采样可以被看成是Metropolis-Hastings算法的一个特例。而这个特殊之处就在于，使用Gibbs采样时，通常是对多变量分布进行采样的。比如说，现在有一个分布p（z₁，z₂，z₃），可见它是定义在三个变量上的分布。这种分布在实际中还有很多，比如一维的正态分布也是关于其均值和方差这两个变量的分布。在算法的第t步，选出了三个值，，。这时，既然和是已知的，那么可以由此获得一个新的，并用这个新的值替代原始的，而新的可以由下面这个条件分布来获得

接下来同样的道理再用一个新的来替代原始的，而这个新的可以由下面这个条件分布来获得

可见刚刚获得的新值已经被直接用上了。同理，最后用来更新，而新的值由下面这个条件分布来获得

按照这个顺序如此往复即可。下面给出更为完整的吉布斯采样的算法描述。

当然如果要从理论上证明Gibbs采样确实可以得到预期的分布，就应该考察它是否满足细致平衡。但是前面也讲过，Gibbs采样是Metropolis-Hastings算法的一个特例。所以其实甚至无需费力去考察其是否满足细致平衡，如果能够证明吉布斯采样就是Metropolis-Hastings算法，而Metropolis-Hastings算法是满足细致平衡，其实也就得到了我们想要的。下面是证明的过程，其中x_-k表示一个向量（x₁，…，x_k-1，x_k+1，…，x_D），也就是从中剔除了x_k元素。

令x=x₁，…，x_D。

当采样第k个元素，

当采样第k个元素，

以Gibbs采样为基础实现的MCMC在自然语言处理中的LDA里有重要应用。如果读者正在或者后续准备研究LDA的话，这些内容将是非常重要的基础。