马尔可夫链蒙特卡罗方法概述：概率分布采样的经典方法

2019年3月23日

统计学和机器学习里经常要处理这样的概率分布：没有简单的解析形式，或者难以直接采样。马尔可夫链蒙特卡罗（Markov Chain Monte Carlo，MCMC）方法就是为这类问题准备的。本文讲 MCMC 的基本概念、主要算法，以及它在贝叶斯推断中的应用。

MCMC 的基本概念和理论基础

要理解 MCMC，我们首先需要了解其名称中包含的两个关键概念：马尔可夫链和蒙特卡罗方法。

马尔可夫链

马尔可夫链是一种特殊的随机过程，其核心特征是“无记忆性”。具体来说，在一个马尔可夫链中，系统下一个状态的概率分布只依赖于当前状态，而与之前的状态历史无关。这个性质可以用数学语言表示为：

P(X_{n+1} = x | X_n = x_n, X_{n-1} = x_{n-1}, ..., X_1 = x_1) = P(X_{n+1} = x | X_n = x_n)

这里， $X_n$ 表示在时间 $n$ 时系统的状态。

举个简单的例子，想象一个醉汉在一条直线上随机行走。每一步，他都有相同的概率向左或向右移动一个单位距离。他下一步的位置只取决于他当前的位置，而不依赖于他是如何到达当前位置的。这就是一个典型的马尔可夫链。

这个醉汉行走的例子包含了马尔可夫链的几个重要特征：

状态空间：在这个例子中，状态空间是醉汉可能所处的所有位置，可以用整数来表示（例如，0表示起点，正整数表示向右移动的步数，负整数表示向左移动的步数）。状态空间可以是离散的或连续的。
时间离散性：醉汉的移动是一步一步进行的，时间是离散的。每一步对应一个时间点。马尔可夫链可以是时间离散的或时间连续的。
无记忆性（马尔可夫性）：这是马尔可夫链最关键的特征。醉汉下一步的位置只依赖于他当前的位置，而不依赖于他之前的移动历史。无论醉汉是如何到达当前位置的，他下一步的行为都是一样的。
状态转移概率：从一个状态到另一个状态的概率是固定的。在这个例子中，无论醉汉在哪个位置，他向左或向右移动的概率都是相等的（各为 $1/2$ ）。这些概率可以用状态转移矩阵来表示。
平稳分布：在某些条件下，当马尔可夫链运行足够长的时间后，它会收敛到一个平稳分布。在醉汉行走的例子中，如果我们观察足够长的时间，我们会发现醉汉在每个位置上出现的概率会趋于一个固定的分布。
遍历性：如果一个马尔可夫链是遍历的，那么无论从哪个初始状态开始，它最终都会到达任何可能的状态，并且收敛到相同的平稳分布。在醉汉行走的例子中，虽然每一步醉汉只能移动一个单位距离，但经过足够长的时间，他可能会到达直线上的任何位置。

在 MCMC 中，我们构造的马尔可夫链的状态空间就是我们感兴趣的参数空间，而状态转移规则则是由具体的 MCMC 算法（如 Metropolis-Hastings 算法或 Gibbs 采样）来定义的。MCMC 方法的核心思想就是构造一个马尔可夫链，使其平稳分布恰好是我们想要采样的目标分布。

蒙特卡罗方法

蒙特卡罗方法是一类依赖随机采样的数值计算方法。这个方法的核心思想是：通过大量随机样本的统计特性来近似解决问题。

例如，我们可以用蒙特卡罗方法来估计圆周率 $\pi$ 。考虑一个边长为 $2$ 的正方形，其内切一个直径为 $2$ 的圆。我们随机在这个正方形内撒点，然后计算落在圆内的点的比例。当点数足够多时，这个比例会接近于 $\pi/4$ 。

import random

def estimate_pi(num_points):
    inside_circle = 0
    for _ in range(num_points):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        if x*x + y*y <= 1:
            inside_circle += 1
    return 4 * inside_circle / num_points

print(estimate_pi(10000000))

蒙特卡洛方法也有一些重要性质：

收敛性：根据大数定律，随着样本数量的增加，蒙特卡罗方法的估计会收敛到真实值。
误差估计：蒙特卡罗方法的标准误差通常与 $1/\sqrt{N}$ 成正比，其中 $N$ 是样本数量。这意味着要将误差减半，我们需要增加 $4$ 倍的样本数。
维度无关性：蒙特卡罗方法的收敛速度与问题的维度无关，这使得它在处理高维问题时特别有优势。
并行化：大多数蒙特卡罗算法都可以很容易地并行化，这在现代计算环境中是一个重要的优势。

这里说“蒙特卡洛方法的收敛速度与问题的维度无关”是一种简化的说法。为了达到给定的精度，所需的样本数量不会随着问题维度的增加而呈指数增长。但在实践中，高维问题通常更复杂，可能需要更多的样本来充分探索整个空间。另外，在某些高维问题中，虽然收敛速度在理论上不变，但估计的方差可能随维度增加而增大，这实际上会影响到所需的样本数量。
而之所以蒙特卡洛方法很容易并行化主要有以下原因:
独立采样：蒙特卡罗方法的核心是独立的随机采样。每个样本的生成通常不依赖于其他样本，这种独立性使得我们可以同时在多个处理器或计算节点上生成样本。
可加性：蒙特卡罗估计通常涉及对大量独立样本的某种形式的平均或汇总。这种操作具有可加性，意味着我们可以将总体任务分解为多个子任务，然后将结果合并。

在 MCMC 方法中，蒙特卡罗的思想被用来从复杂的概率分布中进行采样。通过构造一个特殊的马尔可夫链，我们可以生成服从目标分布的样本，然后使用这些样本来估计各种统计量或进行积分计算。

MCMC 的核心思想

MCMC 方法的核心思想是通过构造一个特定的马尔可夫链来从复杂的概率分布中进行采样。它结合了马尔可夫过程的状态转移机制和蒙特卡罗方法的随机采样原理。

这个思想可以拆成六步：

问题定义：想象我们面对一个复杂的概率分布。这个分布可能是高维的，可能没有简单的数学表达式，甚至可能只能通过一个复杂的函数按比例计算（即我们知道 $f(x)$ ，但只知道 $p(x) ∝ f(x)$ ，而不知道归一化常数）。我们的目标是从这个分布中获取样本。
马尔可夫链的引入：MCMC 的关键在于构造一个特殊的马尔可夫链。这个链的状态空间与我们的目标分布的支撑集相同。例如，如果我们的目标是一个定义在实数上的分布，那么这个马尔可夫链的状态也是实数。
设计转移规则：构造马尔可夫链的核心是设计其转移规则。这个规则决定了链如何从一个状态移动到下一个状态。MCMC 的关键在于，它设计了一种特殊的转移规则，使得当链运行足够长时间后，其状态分布会收敛到我们的目标分布。
收敛到目标分布：MCMC 方法的核心在于，无论我们从何种初始状态开始，经过足够长的时间后，马尔可夫链的状态分布都会收敛到目标分布。这就是所谓的“遍历性”。这个过程可以这样理解：起初，链可能在状态空间中任意游走，但随着时间推移，它会逐渐遗忘初始状态，开始根据目标分布的形状在重要区域间跳转。
从平稳分布中采样：当马尔可夫链达到平稳状态后（通常称为“燃烧期 burn-in period”之后），我们就可以将链上的状态视为来自目标分布的样本。这些样本可能会有一定的自相关性，但它们的分布会近似于目标分布。
蒙特卡罗估计的应用：有了这些样本，我们就可以应用蒙特卡罗方法来估计目标分布的各种统计量。例如，我们可以用样本平均来估计分布的期望，用样本方差来估计分布的方差等。

MCMC 的适用条件只有一个：只要能按比例计算目标分布的概率密度（或质量），就能构造一个 MCMC 算法从中采样，不管这个分布本身多复杂。复杂贝叶斯模型和高维统计推断都落在这个条件里。

Metropolis-Hastings 算法详解

Metropolis-Hastings（MH）算法是最基本也是最重要的 MCMC 算法之一。它提供了一种通用的方法来构造收敛到指定平稳分布的马尔可夫链。

算法步骤

$\pi(x)$ 是我们想要采样的目标分布（Target Distribution），我们通常只能计算这个分布的非标准化形式。提议分布（Proposal Distribution）是我们用来生成新的候选样本的分布，通常记为 $q(x'|x)$ 。提议分布应该易于采样，常见的选择包括高斯分布或均匀分布。因为直接从复杂的目标分布采样通常是困难或不可能的，所以我们需要提议分布来探索参数空间。

MH 算法的步骤如下：

选择一个初始状态 $x_0$ 。
对于 $t = 0, 1, 2, ...$ $t = 0, 1, 2, ...$ :
1. 从一个提议分布 $q(x'|x_t)$ 中生成一个候选状态 $x'$ 。
2. 计算接受概率： $\alpha = min(1, \frac{\pi(x') \cdot q(x_t|x')}{\pi(x_t) \cdot q(x'|x_t)})$
3. 生成一个 $0$ 到 $1$ 之间的随机数 $\mu$ 。
4. 如果 $\mu < \alpha$ ，接受这个移动，令 $x_{t+1} = x'$ ；否则，保持原状态， $x_{t+1} = x_t$ 。

算法原理解析

MH 算法的关键是接受概率 $\alpha$ 。这个概率平衡了目标分布的相对概率密度和提议分布的不对称性。

目标分布的相对概率密度（ $\pi(x') / \pi(x_t)$ ）引导算法向高概率区域移动，确保我们更频繁地采样到目标分布的重要部分。
提议分布的不对称性校正（ $q(x_t|x') / q(x'|x_t)$ ）消除了由于提议分布选择而可能引入的偏差，使得算法的结果不依赖于特定的提议分布。

这两个因素的平衡确保了细致平衡条件的满足，这是马尔可夫链收敛到目标分布 $\pi(x)$ 的充分条件。没有这种平衡，算法会产生有偏的样本，无法正确表示目标分布。

细致平衡条件（Detailed Balance Condition）可以表述为：对于任意状态 $x$ 和 $y$ ，
$\pi(x) \times P(x → y) = \pi(y) \times P(y → x)$
这里， $\pi(x)$ 和 $\pi(y)$ 分别表示系统在状态 $x$ 和状态 $y$ 的平稳分布概率， $P(x → y)$ 是从状态 $x$ 转移到状态 $y$ 的概率。这个等式表明，在平衡状态下，系统从状态 $x$ 转移到状态 $y$ 的总体流量与从状态 $y$ 转移到状态 $x$ 的总体流量相等。这确保了无论初始状态如何，生成的样本序列最终将代表目标分布。

Gibbs 采样

Gibbs 采样是另一种重要的 MCMC 方法，特别适用于多维问题，可以看作是 MH 算法的一个特例，但其独特的采样策略使它在某些情况下特别有效。

算法步骤

Gibbs 采样的核心思想是：每次只更新多维分布中的一个维度，而保持其他维度不变。这种方法利用了条件分布通常比联合分布更容易采样的特点。

假设我们要从 $n$ 维联合分布 $\pi(x_1, ..., x_n)$ 中采样。Gibbs 采样的步骤如下：

选择一个初始状态 $x^{(0)}=(x_1^{(0)}, ..., x_n^{(0)})$ 。
对于 $t = 1, 2, ...$ :
对于 $i = 1$ 到 $n$ :
1. 从条件分布 $\pi(x_i | x_1^{(t)}, ..., x_{i-1}^{(t)}, x_{i+1}^{(t-1)}, ..., x_n^{(t-1)})$ 中采样 $x_i^{(t)}$ 。
2. 更新状态： $x^{(t)}=(x_1^{(t)}, ..., x_{i}^{(t)}, x_{i+1}^{(t-1)}, ..., x_n^{(t-1)})$ 。

我们举一个具体的例子来说明。假设我们要从一个二维正态分布中采样，其中 $X_1$ 和 $X_2$ 有一定的相关性。

初始状态： $x^{(0)}=(x_1^{(0)},x_2^{(0)})=(0,0)$

第一次迭代 ( $t = 1$ ):

对 $X₁$ $X_{1}$ 进行采样：
1. 从 $\pi(x_1 | x_2^{(0)})$ 中采样 $x_1^{(1)}$ 。假设得到 $x_1^{(1)}=0.5$
2. 更新状态： $x^{(1)} = (0.5, 0)$
对 $X_2$ $X_{2}$ 进行采样：
1. 从 $\pi(x_2 | x_1^{(1)})$ 中采样 $x_2^{(1)}$ 。假设得到 $x_2^{(1)} = -0.3$
2. 更新状态： $x^{(1)} = (0.5, -0.3)$

第二次迭代 ( $t = 2$ ):

对 $X₁$ $X_{1}$ 进行采样：
1. 从 $\pi(x_1 | x_2^{(1)})$ 中采样 $x_1^{(2)}$ 。假设得到 $x_1^{(2)}=0.2$
2. 更新状态： $x^{(2)} = (0.2, -0.3)$
对 $X_2$ $X_{2}$ 进行采样
1. 从 $\pi(x_2 | x_1^{(2)})$ 中采样 $x_2^{(2)}$ 。假设得到 $x_2^{(2)} = 0.1$
2. 更新状态： $x^{(2)} = (0.2, 0.1)$

这个过程不断重复，每次迭代都更新所有维度。随着迭代次数增加，生成的样本将逐渐收敛到目标二维正态分布。

关键点是，每次我们只更新一个维度，而这个维度的新值是基于其他维度的当前值采样得到的。这就是 Gibbs 采样的核心思想：通过不断地从条件分布中采样来逼近联合分布。

算法原理解析

Gibbs 采样自动满足细致平衡条件，不需要显式的接受或拒绝步骤。回顾一下细致平衡条件： $\pi(x) \cdot P(x → y) = \pi(y) \cdot P(y → x)$ ，在 Gibbs 采样中，每一步的移动都是沿着条件分布进行的。这个条件分布不是任意选择的，而是目标分布在当前维度上的精确表示。当我们在第 $i$ 个维度上采样时，我们使用的条件分布是 $\pi(x_i | x_1, ..., x_{i-1}, x_{i+1}, ..., x_n)$ 。这个条件分布实际上是从完整的联合分布 $\pi(x_1, ..., x_n)$ 中导出的。因此，当我们从这个条件分布中采样时，我们实际上是在目标分布的一个“切片”上移动。因为我们总是沿着目标分布的精确条件分布移动，所以每一步都自动满足细致平衡条件。我们不需要额外的接受或者拒绝步骤来纠正偏差。

Gibbs 采样利用条件分布的特性，把复杂的多维采样问题转化成一系列一维采样问题。贝叶斯推断和图模型是它最常见的应用场景。

MCMC 在贝叶斯推断中的应用

MCMC 方法在贝叶斯统计中扮演着关键角色，使得复杂的后验分布推断成为可能。

贝叶斯推断简介

贝叶斯推断的核心思想是将参数视为随机变量，并通过观测数据来更新我们对参数的信念。这个过程可以用贝叶斯定理来描述：

P(\theta|D) = \frac{(P(D|\theta) * P(\theta))}{P(D)}

这里， $θ$ 是模型参数， $D$ 是观测数据， $P(\theta)$ 是参数的先验分布，表示我们在看到数据之前对参数的信念， $P(D|\theta)$ 是似然函数，描述给定参数下观测到数据的概率， $P(\theta|D)$ 是后验分布，表示我们在看到数据后对参数的更新信念， $P(D)$ 是证据，通常作为归一化常数。

MCMC在贝叶斯推断中的作用

在许多实际问题中，后验分布 $P(θ|D)$ 往往没有解析形式，或者难以直接计算。这就是 MCMC 方法发挥作用的地方。MCMC 允许我们：

直接从后验分布中采样，而不需要显式计算后验分布。
估计后验分布的各种统计量，如均值、方差、分位数等。
进行后验预测。

下一步: MCMC 的自我"链"接

标准 MCMC 有一个未解决的问题：提议分布是固定的。在高维空间或结构复杂的分布中，固定的提议分布会导致收敛缓慢或采样效率低下，而合适的尺度又很难事先猜到。

自适应 MCMC 的思路是在运行过程中利用已经获得的样本信息来改进采样策略。下一篇介绍它的基本原理和发展历程、自适应 Metropolis 算法（AM）的工作机制，以及自适应提议分布如何提高采样效率。

参考文献

Robert, C. P., & Casella, G. (2004). Monte Carlo statistical methods (Vol. 2). New York: Springer.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian data analysis. CRC press.
Brooks, S., Gelman, A., Jones, G., & Meng, X. L. (Eds.). (2011). Handbook of markov chain monte carlo. CRC press.

差分进化算法简介：优化算法的创新自适应 MCMC 方法：提高采样效率的尝试