马尔可夫链蒙特卡罗(MCMC)方法是贝叶斯推理中常用的一种非常强大的蒙特卡罗方法。

“经典”蒙特卡罗方法依赖于由独立观测数据组成的计算机生成样本，而MCMC方法则用于生成相关观测数据序列。因为这些序列是马尔可夫链，这也是该方法名称的由来。

蒙特卡罗

要理解MCMC，我们需要熟悉蒙特卡罗方法的基础知识。

当我们不能用解析的方法计算出来时，我们使用蒙特卡罗方法来近似一个随机变量X的概率分布的特征(例如，它的期望值)。

通过计算机，我们从X的分布中生成独立抽取的样本$x_1, \ldots, x_T$。

然后，我们使用样本的经验分布（ empirical distribution）来进行我们的计算。经验分布是离散概率分布，将概率$1/T$分配给每个值$x_1, \ldots, x_T$。

例子 如果我们需要计算X的期望值，我们用经验分布的期望值近似它，这就是$x_1, \ldots, x_T$的样本均值。

例子 如果我们需要计算X低于某个阈值的概率，我们可以用计算机生成的样本中低于阈值的观测值的百分比来近似它。

例子 为了近似X的方差，我们计算$x_1, \ldots, x_T$的样本方差。

plug-in的原则

一般来说，无论我们想对X的真实分布进行什么计算，都要近似于对经验分布进行的相应计算。后者很容易在计算机上执行，因为它涉及一个具有有限支持的离散分布。

这种方法通常工作得很好，在某种意义上，随着样本量$T$的增加,近似收敛于真实量。这就是所谓的plug-in原则。

马尔可夫链基础

MCMC方法的工作原理类似于标准的蒙特卡罗方法，但有一点不同:计算机生成的$x_1, \ldots, x_T$不是独立的，但它们是序列相关的。

具体来说，它们是$T$随机变量$X_1$，…， $X_{T}$组成一个马尔可夫链。

马尔可夫性质

随机序列$x_1, \ldots, x_T$是马尔可夫链，当且仅当给定当前值$X_{t}$，对于任何正整数$k$和$n$，未来的观察值$X_{t+n}$与过去的值$X_{t-k}$有条件独立。

这个性质被称为马尔可夫性质，它表示这个过程是无记忆的:链的未来值的概率分布只依赖于它的当前值$X_{t}$，而不管该值是如何达到的(即不管链在过去遵循的路径是什么)。

条件分布和无条件分布

由于马尔可夫性质，我们基本上只需要两个东西来分析一条链:

$X_{t+n}= x_{t+n} $的条件概率(或密度)，假设$X_{t}= x_{t} $，表示为

$$ f\left(x_{t+n} \mid x_t\right) $$

$X_{t}= x_{t} $的无条件概率，用 $$ f\left(x_t\right) $$

渐近独立(Asymptotic independence)

虽然在一般的马尔可夫链中并不是如此，但MCMC方法生成的链具有以下性质:

这个性质意味着当n变大时$f\left(x_{t+n} \mid x_t\right)$收敛到$f\left(x_t\right)$。

这一点非常重要。我们通常以某种任意的方式生成链的第一个值($x_{1}$)，但我们知道，当$t$变大时，链项的分布受初始选择的影响越来越小。换句话说，无论我们如何选择$x_{1}$， $f\left(x_{t+n} \mid x_t\right)$都收敛于相同的分布$f\left(x_t\right)$。

目标分布(Target distribution)

在MCMC方法生成的链中，不仅$f\left(x_{t+n} \mid x_t\right)$收敛于$f\left(x_t\right)$，而且，随着$t$的变化，分布$f\left(x_t\right)$彼此之间变得几乎相同。它们收敛于所谓的链的平稳分布。

而且，平稳分布与目标分布重合，也就是我们要从中提取样本的分布。

也就是说，$t$越大，$f\left(x_{t+n} \mid x_t\right)$和$f\left(x_t\right)$越接近目标分布。

黑箱方法(A black-box approach)

让我们总结一下到目前为止所介绍的内容。

我们可以把MCMC算法想象成一个接受两个输入的黑盒子:

  1. 初始值$x_{1}$;

  2. 目标分布。

输出是一系列值$x_1, \ldots, x_T$。

虽然序列的第一个值是从可能与目标分布非常不同的分布中提取的，但这些分布与目标之间的相似性逐渐增加:当$t$很大时，$x_{t}$是从与目标非常相似的分布中提取的。

老化样本(Burn-in sample)

由于目标分布和链的第一项分布之间的差异，通常的做法是丢弃MCMC样本的第一次抽取。

例如，如果我们有110,000张图，我们丢弃前10,000张，保留剩下的100,000张。

被丢弃的一组称为老化岩本。

通过丢弃老化样品，我们消除了分布远离目标分布的样本，而保留了分布接近目标分布的样本。

这减少了用MCMC样本执行的任何蒙特卡罗近似的偏差。

相关性和有效样本量

丢弃老化样本后，我们从分布中抽取的样本与目标非常相似。

然而，我们仍然有一个问题:我们的绘图不是独立的。

缺乏独立性意味着什么?

要理解它们，我们需要回到标准的蒙特卡洛方法。

标准蒙特卡罗模拟的准确性取决于样本量$T$: $T$越大，近似值越好。

在MCMC模拟的情况下，我们需要使用有效样本量的概念:$T$相关的观测值等价于更少的独立观测值。

例如，1000个相互依赖的观测值可以等价于100个独立的观测值。在这种情况下，我们说有效样本容量等于100。

粗略地说，相邻观测值之间的相关性越高，有效样本量越小，MCMC近似越不准确。

这就是为什么在MCMC模拟中，大部分的努力都致力于尽可能地减少相关性，以达到令人满意的有效样本量。

Metropolis-Hastings算法

最流行的MCMC算法之一是Metropolis-Hastings (M-H)算法。

用$f(x)$表示目标分布的密度(或质量)函数，我们希望从中提取绘制序列的分布。

例如，$f$可以是贝叶斯推理中的后验密度。

用$q\left(x \mid x^{\prime}\right)$表示一个我们能够从中提取计算机生成样本的条件分布(例如，一个均值为$x^{prime}$的多元正态分布)。

生成的样本$x$、我们所条件的值$x^{prime}$和目标分布的参数$f$都具有相同的维数。

算法

M-H算法从属于目标分布支持的任意值$x_{1}$开始，生成后续值$x_{t}$如下所示:

收敛

在满足一定技术条件的情况下，可以证明$\{x_t\}$的分布收敛于目标分布$f(x)$

此外, $$ \frac{1}{T} \sum_{t=1}^T g\left(x_t\right) 近似趋近于 \int g(x) f(x) d x $$ 其中$g$是期望值$\int|g(x)| f(x) d x$存在且有限的任何函数。

换句话说，一个强大的大数定律(一个遍历定理)适用于样本均值$\frac{1}{T} \sum_{t-1}^T g\left(x_t\right)$。

M-H算法的主要优点

Metropolis-Hastings算法的强大之处在于你只需要知道函数$f$和它乘上一个乘法常数即可。

例如，在贝叶斯推断中，非常常见的是知道一个乘法常数的后验分布，因为似然和先验已知，但边际分布不知道。

假设 $f(x)=\operatorname{ch}(x)$ 我们知道$h(x)$，但不知道$c$。

M-H算法中的接受概率为 $$ \begin{aligned} p_t &=\min \left(\frac{f\left(y_t\right)}{f\left(x_{t-1}\right)} \frac{q\left(x_{t-1} \mid y_t\right)}{q\left(y_t \mid x_{t-1}\right)}, 1\right) \\ &=\min \left(\frac{\operatorname{ch}\left(y_t\right)}{\operatorname{ch}\left(x_{t-1}\right)} \frac{q\left(x_{t-1} \mid y_t\right)}{q\left(y_t \mid x_{t-1}\right)}, 1\right) \\ &=\min \left(\frac{h\left(y_t\right)}{h\left(x_{t-1}\right)} \frac{q\left(x_{t-1} \mid y_t\right)}{q\left(y_t \mid x_{t-1}\right)}, 1\right) \end{aligned} $$ 因此，接受概率是唯一依赖于$f$的量，可以在不知道常数$c$的情况下计算出来。

这就是Metropolis-Hastings算法的美妙之处:即使我们不完全知道该分布的密度，我们也可以从该分布生成绘图!

吉布斯抽样(Gibbs sampling)

另一种流行的MCMC算法是所谓的吉布斯采样算法。

充分条件

假设我们想要生成一个随机向量$x_{ullet}$的采样，它具有联合概率密度 $$ f\left(x_{\bullet}\right)=f\left(x_{\bullet}, 1, \ldots, x_{\bullet, B}\right) $$

给定块$x_{\cdot,b}$，用$x_{\cdot，-b}$表示除$x_{\cdot,b}$外的所有块组成的向量。

假设我们能够从$B$条件分布$x_{\bullet}, 1, \ldots, x_{\bullet, B}$中提取条件分布 $$f\left(x_{\bullet, 1} \mid x_{\bullet,-1}\right), \ldots f\left(x_{\bullet . B} \mid x_{\bullet,-B}\right)$$

在MCMC术语中，这些被称为完全条件分布。

算法

吉布斯采样算法从任意向量开始 $$x_1=\left[\begin{array}{lll} x_{1,1} & \ldots & x_{1, B} \end{array}\right]$$

并生成后续值$x_{t}$如下所示:

注意，在迭代$t$时，当您提取$b$-th块时，块1到$b-1$的最新提取是在迭代$t$中已经提取的，而块$b+1$到$b$的最新提取是在前一个迭代($t-1$)中提取的。

收敛

可以证明吉布斯抽样是Metropolis-Hastings抽样的一个特例。因此，就像在M-H中一样，抽样的分布收敛于目标分布$f(x)$。此外，样本均值的一个遍历定理$$ \frac{1}{n} \sum_{i=1}^n g\left(x_i\right) $$成立。

参考资料：