Markov Chain Monte Carlo (MCMC)

约 2601 个字 预计阅读时间 9 分钟 总阅读量 次

顺带帮你复习贝叶斯公式！

MCMC（Markov Chain Monte Carlo，马尔可夫链蒙特卡洛方法）是现代统计学、物理学和机器学习中非常核心的一类算法。是一种==在无法直接计算积分的情况下==，通过构建马尔可夫链来从复杂概率分布中进行采样，从而近似计算期望值或积分的方法。

What is it?

拆解开来最好理解 MCMC：

1. 蒙特卡洛（Monte Carlo）：
   • 指使用随机数（采样）来解决计算问题的方法。如想求一个复杂形状的面积（积分），你可以向这个形状里随机扔很多石子，通过统计石子落在形状内的比例来估算面积。
   • 在统计学中，我们通常用它来计算期望值，其中 \(x_i\) 是从分布 \(p(x)\) 中抽取的样本。：

\[E[f(x)] \approx \frac{1}{N} \sum_{i=1}^{N} f(x_i)\]

1. 马尔可夫链（Markov Chain）：
   • 指一个随机过程，在这个过程中，下一个状态只取决于当前状态，而与过去的状态无关（“无记忆性”）。
   • 关键性质：如果在合适的条件下运行足够长的时间，马尔可夫链会达到一个平稳分布（Stationary Distribution）。即不管初始状态在哪，最终状态的分布会收敛到一个特定的分布。

MCMC 的结合点

我们想要从一个非常复杂的分布 \(p(x)\) 中采样，但直接采样太难了。于是，我们巧妙地设计一个马尔可夫链，使得这个链的平稳分布恰好就是我们需要采样的目标分布 \(p(x)\)。然后，我们让这个链运行一段时间，记录下它经过的状态，这些状态就是我们想要的样本。

---

2. MCMC 为了解决什么样的问题？

MCMC 主要解决的是高维积分和归一化常数无法计算的问题，这在 贝叶斯推断（Bayesian Inference） 中最为常见。

贝叶斯推断

贝叶斯推断与贝叶斯概率密切相关，其通过处理先验、证据和似然来计算后验概率。给定某个事件B，事件A发生的概率是多少？这可以通过贝叶斯著名的公式来回答：\(P(A/B)=\dfrac{P(B/A)P(A)}{P(B)}\)。

• \(P(A/B)\) 是 后验概率(Posterior)。这是我们希望计算的。
• \(P(B/A)\) 是 似然 (Likelihood)。假设A已发生，B发生的可能性有多大。
• \(P(A)\) 是 先验概率 (prior)。在不考虑证据的情况下，事件\(A\)本身的可能性。
• \(P(B)\) 是 证据(Evidence)。在不考虑事件的情况下，证据\(B\)本身的可能性。

我们主要关注贝叶斯公式的以下具体形式：

\(P(\theta/D)=\dfrac{P(D/\theta)P(\theta)}{P(D)}\)，其中\(P(\theta/D)\) 是 后验概率，\(P(D/\theta)\) 是 似然，\(P(\theta)\) 是 先验概率，而\(P(D)\) 是 证据。

我们的目标是找到最可能的参数\(\theta\)的分布，这些参数能够解释观测数据D。

计算这些概率中的某些部分可能很繁琐，尤其是证据\(P(D)\)。此外，还可能遇到其他问题，例如确保共轭性（本文不深入探讨此问题）。幸运的是，一些技术，即MCMC，允许我们从后验分布中抽样，并对参数进行分布推断，而无需担心计算证据或共轭性问题。

问题出在分母 \(P(X)\) 上：

\[P(D) = \int P(D | \theta) P(\theta) \, d\theta\]

在参数 \(\theta\) 维度很高时（例如神经网络中有数百万个参数，或者复杂的统计模型），这个积分是根本无法算出解析解的，数值积分（如梯形法则）在维数灾难下也会失效。

MCMC 的解决方案：

MCMC 允许我们在不知道分母 \(P(X)\) 的情况下，从分子 \(P(X | \theta) P(\theta)\)（即未归一化的后验概率）中采样。只要能比较两个点的相对概率大小，MCMC 就能工作。

换句话说，MCMC 允许我们从任何无法直接抽样的分布中进行抽样。它可用于从参数的后验分布中抽样。该方法已在许多应用中取得了巨大成功，例如在给定一组观测数据和先验信息的情况下计算参数的分布，以及在物理学和数字通信中计算高维积分。

总而言之：在给定一组观测数据和一个先验信息的情况下，它可以用于计算参数的分布。

---

3. 基本流程

为何可以通过“比较相对概率大小”来绕过“无法直接计算绝对概率”的难题。

主要用到以下的内容，用符号表示为：

• \(\theta\) (Theta): 我们想要推断的参数（例如：某种统计模型的参数，或者神经网络的权重）。
• \(D\) (Data): 我们观测到的数据。
• \(P(\Theta=\theta | D)\): 后验分布（Posterior）。这是我们的终极目标，即在看到数据 \(D\) 后，参数 \(\theta\) 的概率分布。这通常是 "intractable"（难以处理的），因为分母无法计算。
• \(Q(\theta' | \theta)\): 提议分布（Proposal Distribution / Transition Model）。
  • 这是一个用来帮助我们在参数空间中“随机游走”的工具。通常选择简单的分布，比如高斯分布 \(\mathcal{N}(\theta, \sigma)\)。
  • 作用：基于当前位置 \(\theta\)，给出一个建议的下一跳位置 \(\theta'\)。
• \(f\): 非归一化的目标函数，每一个新样本的“可能性”由函数 \(f\) 决定。这就是为什么 \(f\) 必须与我们要采样的后验分布成正比。通常选择 \(f\) 为能够表达这种比例关系的概率密度函数。
  • 文本解释：\(f\) 正比于后验分布。通常 \(f = \text{Likelihood} \times \text{Prior}\)（似然 \(\times\) 先验）。

MH 算法最巧妙的地方在于它不需要算出准确的后验概率，只需要比较两个位置的概率大小。

A. 贝叶斯公式的简化

事实上贝叶斯公式可以做一个（相对）扩展：

\[\frac{P(\theta' | D)}{P(\theta | D)} = \frac{\frac{P(D|\theta')P(\theta')}{P(D)}}{\frac{P(D|\theta)P(\theta)}{P(D)}}\]

这里发生了一个“魔法”消除：分母 \(P(D)\)（证据 Evidence）被上下消掉了。 $$ Ratio = \frac{P(D|\theta')P(\theta')}{P(D|\theta)P(\theta)} = \frac{\text{新位置的似然} \times \text{新位置的先验}}{\text{旧位置的似然} \times \text{旧位置的先验}} $$

这意味着我们只需要算出分子（Likelihood \(\times\) Prior），就可以决定是否接受新样本。

---

B. 接受率 \(P(\text{accept})\)

是否跳转的规则：

\[ P(\text{accept}) = \min\left( 1, \frac{f(\text{新})}{f(\text{旧})} \right) \]

• 情况 1：如果新位置 \(\theta'\) 概率更高（比值 \(>1\)）
  • 必定接受 (\(P=1\))。这就像爬山，看到更高的地方，立刻爬上去。
• 情况 2：如果新位置 \(\theta'\) 概率更低（比值 \(<1\)）
  • 以一定概率接受。比如比值是 0.3，我们就有 30% 的机会接受它，70% 的机会拒绝它。
  • 为什么要接受更差的位置？ 为了探索（Exploration）。如果只往高处走，就会陷入局部最优（Local Optima）。接受坏样本允许算法跳出并找到全局最优分布。

---

C. 合起来看

Metropolis-Hastings 使用 \(Q\) 在分布空间中进行随机游走，并根据样本的似然（likelihood）来接受或拒绝跳转到新的位置。这种“无记忆（memoriless）”的随机游走正是 MCMC 中“马尔可夫链（Markov Chain）”的部分。

为了得到参数的新位置，只需取我们当前的位置 \(\theta\)，并提议（propose）一个新的位置 \(\theta^\prime\)，它是从分布 \(Q(\theta^\prime/\theta)\) 中抽取的随机样本。这通常是一个对称分布。例如，一个均值为 \(\theta\) 且具有一定标准差 \(\sigma\) 的正态分布：\(Q(\theta^\prime/\theta) = \mathcal{N}(\theta, \sigma)\)。

为了决定 \(\theta^\prime\) 是被接受还是被拒绝，必须为每一个新的 \(\theta^\prime\) 值计算以下比率： \(\dfrac{P(\theta^\prime/D)}{P(\theta/D)}\)。利用贝叶斯公式，这可以很容易地重新表述为：\(\dfrac{P(D/\theta^\prime)P(\theta^\prime)}{P(D/\theta)P(\theta)}\)（证据项 \(P(D)\) 在除法过程中被直接约掉了）。

\(\dfrac{P(D/\theta^\prime)P(\theta^\prime)}{P(D/\theta)P(\theta)}\) 也等同于 \(\dfrac{\prod_i^nf(d_i/\Theta=\theta^\prime)P(\theta^\prime)}{\prod_i^nf(d_i/\Theta=\theta)P(\theta)}\)。

似然是什么？

这里的 \(\prod_i^nf(d_i/\Theta=\theta^\prime)\) 就是“似然”，表示“现有的数据觉得这个参数有多靠谱”，而 \(P(\theta)\) 则是我们给出的“先验”，我们觉得这个数据有多靠谱；

似然（Likelihood） \(P(D|\theta)\) 不是通过“之前的采样记录”或“历史数据”统计出来的。 它是通过你建立的数学模型函数直接根据当前这一个 \(\theta\) 和 所有的观测数据 \(D\)* 现算出来的。

似然函数本质上是你在回答这个问题：“假设（如果）世界的真实参数就是 \(\theta\)，那么我看到眼前这堆数据 \(D\) 的概率有多大？”

你需要两样东西来计算它： 1. 固定的观测数据 \(D\)：这是你手里拿到的真实数据，比如这里有 100 个硬币的抛掷结果。这些数据在整个 MCMC 过程中是永远不变的。 2. 假设的模型公式：你必须预先根据业务背景定义一个概率分布公式。

---

一般化公式：\(\prod f(d_i | \theta)\) 的理解：

现在回到那个复杂公式： $$ \prod_i^n f(d_i | \Theta=\theta) $$

• \(\prod\) (连乘)：把每一个数据点算出的概率乘起来。
• \(d_i\)：第 \(i\) 个真实的观测数据（比如第 \(i\) 个人的身高）。
• \(\theta\)：当前这一步 MCMC 链条假设的参数值（比如假设平均身高 170）。
• \(f\)：假设的概率密度函数（比如正态分布函数）。

操作步骤： 1. 你手里有 1000 个数据点 \(D\)。 2. MCMC 给你一个参数 \(\theta\)（比如 \(\mu=170, \sigma=10\)）。 3. 你把这 1000 个数据点，逐个代入到正态分布公式 \(PDF(x; 170, 10)\) 中，算出 1000 个概率密度值。 4. 把这 1000 个值乘起来（或者取对数加起来）。 5. 这就是Total Likelihood。

先验的 P 是什么？

Prior（先验）\(P(\theta)\) 不是通过数据“算”出来的，而是由你（主要构建者）在开始之前“选”或者“定义”出来的。

在 MCMC 运行的代码中，它仅仅是一个函数。当你给它一个参数 \(\theta\) 时，它返回一个概率密度值（Density Value）。

• 最常见的情况：均匀先验（Uniform Prior）

如果你对参数 \(\theta\) 一无所知，或者认为它取什么值都一样公平，你会假设它是“均匀分布”的。

• 定义：在某个范围内（比如 -100 到 100），\(P(\theta)\) 是一个常数 \(C\)。
• 计算：无论 \(\theta\) 是什么，函数都返回 \(C\)。
• 在公式中的表现：

当我们计算接受率 Ratio 时：

\[\text{Ratio} = \frac{\text{Likelihood}' \times P(\theta')}{\text{Likelihood} \times P(\theta)} = \frac{\text{Likelihood}' \times C}{\text{Likelihood} \times C} = \frac{\text{Likelihood}'}{\text{Likelihood}}\]

也就是所谓的“先验项通常被约掉”：在这种情况下，你确实不需要计算先验，因为它上下抵消了。

• 带有信息的情况：正态先验（Normal Prior）：当你认为参数大概率在某个值附近时

假设你在估算一个人的身高 \(\theta\)。虽然你还没看数据，但根据常识，你觉得身高大概率在 170cm 左右，不太可能是 3米。既然如此，你可以定义先验为均值 170 的正态分布。

• 定义：\(P(\theta) = \mathcal{N}(170, 10)\)。
• 如何计算：这里就用到了正态分布的概率密度函数（PDF）公式。 $$ P(\theta) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{{-\frac{(\theta-\mu)}2}{2\sigma^2}} $$ 在写代码时，你不需要手写这个复杂公式，通常直接调用库函数。
  • 比如（Python伪代码）：prior_score = scipy.stats.norm(170, 10).pdf(theta)
• 在 MCMC 中的作用： 如果 MCMC 游走到 \(\theta = 300\)（3米高）：
  • \(P(300)\) 会算出一个极小的值（接近0）。
  • 这意味着 \(P(\theta')\) 非常小，会导致整个 \(\text{Ratio}\) 很小。
  • 结果：算法会倾向于拒绝这个“不靠谱”的参数值。

于是：

\[\begin{equation} P(\text{接受}) = \begin{cases}\dfrac{\prod_i^nf(d_i/\Theta=\theta^\prime)P(\theta^\prime)}{\prod_i^nf(d_i/\Theta=\theta)P(\theta)}, & \text{若新位置概率} < \text{旧位置概率} \\  1, & \text{若新位置概率} \ge \text{旧位置概率} \end{cases} \end{equation}\]

这意味着，如果 \(\theta^\prime\) 比当前的 \(\theta\) 更有可能（概率更高），那么我们总是接受 \(\theta^\prime\)。如果它的可能性低于当前的 \(\theta\)，那么我们可能会接受它，也可能会拒绝它，且其被接受的概率随着可能性的降低而减小。

---

D. 整理一下伪代码：

• 指南针 \(Q(\theta'|\theta)\)：即 提议分布。
  • 文中设定为一个对称分布，比如以当前位置为中心的正态分布 \(\mathcal{N}(\theta, \sigma)\)。这就好比你在当前位置，随机向周围迈一步。
• 高度计 \(f\)：即 未归一化的后验概率。
  • \(f(\theta) \propto P(D|\theta)P(\theta)\)。
  • 它由两部分组成：似然（\(\prod f(d_i|\theta)\)，即数据觉得这个参数有多靠谱）乘以 先验（\(P(\theta)\)，即我们要觉得这个参数有多靠谱）。

假设我们现在处于第 \(t\) 步，当前参数位置是 \(\theta\)。

Step 1: 提议（The Propose Step） * 动作：基于当前位置 \(\theta\)，利用 \(Q\) 生成一个候选位置 \(\theta'\)。 * 原理：这就是文本中提到的“无记忆随机游走（Memoriless Random Walk）”。 * 例子：\(\theta' = \theta + \text{随机噪声}\)。

Step 2: 评估比率（The Ratio Step） * 动作：计算“新位置”和“旧位置”的概率比值。 * 公式： $$ \text{Ratio} = \frac{\text{新位置的得分}}{\text{旧位置的得分}} = \frac{Likelihood(\theta') \times Prior(\theta')}{Likelihood(\theta) \times Prior(\theta)} $$ * 关键点：文本特意强调，\(P(D)\) 被约掉了，所以我们不需要算出具体的概率值，只需要算出相对倍数。

Step 3: 决策（The Accept/Reject Step） 这是文本中最核心的逻辑，分为两种情况：

• 情况 A：新位置更好 (\(\text{Ratio} \ge 1\))

  • 判定：\(P(\text{accept}) = 1\)。
  • 结果：必须接受。直接跳到 \(\theta'\)。
  • 直觉：既然新参数比旧参数更能解释数据，当然要留下来。

• 情况 B：新位置更差 (\(\text{Ratio} < 1\))

  • 判定：\(P(\text{accept}) = \text{Ratio}\)。
  • 结果：以概率接受。我们抛一枚不均匀的硬币（在这个 [0,1] 区间内生成一个随机数 \(r\)）。
    • 如果 \(r < \text{Ratio}\)：虽然它变差了，但还是接受。
    • 如果 \(r \ge \text{Ratio}\)：拒绝，原地不动（保持 \(\theta\) 不变，并将 \(\theta\) 再次记入样本序列）。
  • 直觉：虽然新位置不好，但如果它只比现在差一点点（Ratio=0.9），我们要有很大几率过去；如果它差得离谱（Ratio=0.001），我们过去的几率就微乎其微。这正是为了防止陷入局部最优，保持探索性。

通过不断重复 游走 -> 比较 -> 依概率接受，这个链条产生的 \(\theta\) 序列，最终会在概率高的地方出现得密集，在概率低的地方出现得稀疏。这个密度的分布，就完美复现了我们无法直接计算的后验分布 \(P(\theta|D)\)。

---

4. 常见的 MCMC 算法有哪些？

MH 算法是基础，为了提高效率（避免在高维空间中盲目乱走导致拒绝率过高），发展出了更高级的方法：

算法名称	特点	适用场景
Metropolis-Hastings (MH)	最通用的框架，基于拒绝采样。	维度较低，分布形状简单。
Gibbs Sampling	MH 的特例。每次只更新一个维度，固定其他维度。接受率恒为 1（无需拒绝）。	当条件概率 \(p(x_i \| x_{-i})\) 容易写出解析式时（常见于分层贝叶斯模型）。
Hamiltonian Monte Carlo (HMC)	借用物理学中的“哈密顿动力学”。引入动量变量，利用梯度（Gradient）信息来引导采样方向。	高维复杂分布。这是目前概率编程语言（如 PyMC, Stan, Pyro）中的主流算法。
NUTS (No-U-Turn Sampler)	HMC 的自动调参版本。无需手动设置步长和步数。	现代贝叶斯推断的首选默认算法。

---

其他类似的解决方法

MCMC 是一种“采样（Sampling）”方法。如果你的目标是处理贝叶斯推断或复杂分布，还有以下几类替代方案：

A. 变分推断 (Variational Inference, VI)

• 原理：将推断问题转化为优化问题。找一个简单的分布族 \(q(z; \lambda)\)（如高斯分布），通过调整参数 \(\lambda\)，使得 \(q\) 尽可能接近目标后验分布 \(p\)（通常是通过最小化 KL 散度）。
• 对比 MCMC：
  • 速度：VI 通常比 MCMC 快很多，适合大数据集。
  • 精度：VI 是近似解，通常会低估方差；MCMC 理论上只要跑得够久，能得到精确解（无偏）。
  • 应用：VAE（变分自编码器）就是基于 VI；MCMC 常用于需要高精度的科学计算。

B. 拒绝采样 (Rejection Sampling)

• 原理：用一个简单的包络分布覆盖目标分布，采样后决定是否保留。
• 缺点：在高维空间中，接受率会呈指数级下降（几乎所有样本都被拒绝），即“维数灾难”。因此在高维下不可行，必须用 MCMC。

C. 重要性采样 (Importance Sampling)

• 原理：从一个容易采样的分布中采样，然后通过给样本加权重（Weight）来修正偏差。
• 缺点：在高维空间中，权重方差会变得极大，导致只有极少数样本起作用，效率极低。

D. MAP (Maximum A Posteriori)

• 原理：不求整个后验分布，只求后验概率最大的那个点（众数）。等价于加了正则化的最大似然估计。
• 缺点：这是一个点估计，无法给出结果的不确定性（方差/置信区间）。MCMC 能给出完整的分布全貌。

总结

• MCMC 是什么？ 一种构建马尔可夫链来从复杂分布中“游走”并采集样本的方法。
• 核心用途？ 解决贝叶斯推断中分母 \(\int p(x)dx\) 无法计算的难题，从而获得后验分布。
• 流程？ 提议 \(\to\) 比较概率 \(\to\) 接受或拒绝 \(\to\) 循环。
• 进阶？ 高维下首选 HMC/NUTS（利用梯度指引方向）。
• 替代？ 想要速度选 VI（变分推断），想要精确分布选 MCMC。