马尔可夫链蒙特卡洛：藏在现代 AI 底下的 1953 年算法

原文：https://dev.to/lovestaco/markov-chain-monte-carlo-the-1953-algorithm-hiding-under-modern-ai-5cb4（作者 @lovestaco）

有一个算法，1953 年诞生在一台内存比本页 favicon 还小的机器上。

它被用来预报天气。

它被用来拟合黑洞并合的模型。

它潜伏在你导入过的每一个正经贝叶斯统计库的引擎盖底下。

而大多数在职开发者从未听说过它的名字。

它叫马尔可夫链蒙特卡洛（Markov chain Monte Carlo），朋友们都叫它 MCMC。它之所以感觉像个黑盒，是因为人们通常倒着讲它。

他们上来就讲细致平衡、遍历性、平稳分布，等讲到真正巧妙的部分时，你早就把标签页关了。

所以我们顺着讲。

MCMC 是把两个想法拧在一起，而这两个想法都简单到能在酒吧里讲清楚。

想法一：你可以靠扔东西来测量事物

假设我问你圆周率 π 的值，然后收走你的计算器。

你可以把它推导出来。有人推导过，过程很不愉快。

你也可以画一个边长 2R 的正方形，在里面画一个半径 R 的圆，然后蒙着眼睛往正方形上扔飞镖。

一支均匀落在正方形内的飞镖，有一定概率落在圆内，而这个概率正好就是两个面积之比。

圆是 π * R^2。正方形是 4 * R^2。比值是 π / 4。

所以多扔一些飞镖，数一数有多少落在了圆内，乘以四，就得到了 π。

这就是蒙特卡洛方法的全部。它以赌场命名，因为当初在洛斯阿拉莫斯发明它的人正在做中子扩散计算，需要一个代号，而其中一个人的叔叔老是借钱去蒙特卡洛赌钱。

五行代码：

import random

hits = sum(1 for _ in range(10_000_000)
           if random.random()**2 + random.random()**2 <= 1.0)

print(4 * hits / 10_000_000)   # 3.1417...

这段代码里没有人解积分。他们只是数了数。

关键在于——这也是本文剩余部分存在的意义——random.random() 让我们白拿到了独立样本。

我们知道自己在从什么形状里采样。是正方形。从正方形里均匀采样很容易。

而在每一个你真正关心的问题里，你并不知道如何从那个形状里采样。这才是问题所在。

想法二：一个只记得自己此刻在哪里的过程

马尔可夫链是一串状态序列，其中下一个状态只取决于当前状态。

不取决于你是怎么到这里来的。不取决于之前的四十步。只取决于此刻。

这叫作马尔可夫性质（Markov property），简而言之，这条链是"无记忆"的。

教科书例子是天气，那就用教科书例子。三个状态：雨、多云、晴。状态之间的跳转概率是固定的。

从雨出发，有 60% 概率到多云，40% 概率留在原地。从多云出发，有 50% 概率到晴。以此类推。

现在跑一遍：雨、多云、晴、晴、多云。

第五步只参考了第四步，其他什么都没看。它完全不知道这条序列是从雨天开始的。

下面这个性质让它从「只是好玩」变成「真正有用」。

如果这条链能到达每个状态，又不会永久困在任何一个地方，那么它花在每个状态上的时间比例会收敛到一组固定的数字。

跑一千步，你可能有 30% 的时间是雨天。跑一百万步，还是 30%。它稳定下来了。

这组数字就是平稳分布（stationary distribution），而且它不取决于你从哪里出发。

先停下来想一想，因为整个窍门就在这里。

一条马尔可夫链，只要放着不管，就会从某个分布里生成样本。不是某个你事先选好的分布，而是从你随手写下的转移规则里自然冒出来的那个分布。

那么，如果把这个过程倒过来跑呢？

如果你有一个想要的分布，然后你设计转移规则，使它的平稳分布恰好就是这个分布呢？

那你就有了一台机器，可以从一个你从来无法直接采样的分布里，源源不断地吐出样本。

这就是 MCMC。这就是全部想法。剩下的都是工程。

没人能计算的那个分布

先把「一个你无法从中采样的分布」到底是什么说清楚，否则这一切都太空泛了。

假设你在做贝叶斯推断。你有数据，有一个带若干参数的模型，你想知道哪些参数取值与观测结果一致。

非贝叶斯方法会交给你一个最佳拟合数值和一个标准误。

贝叶斯推断交给你的是参数上的一个完整分布，叫后验分布。它既告诉你哪些值是可能的，也告诉你对这些值可以有多大的把握。

这个分布来自贝叶斯定理：四个符号，一个巨大的麻烦。

分子没问题。似然项衡量「这些参数解释我的数据解释得有多好」，这是可以求值的。

先验项回答「看到数据前我相信什么」，这是你自己写下的。

垮掉的是分母。

它叫证据，是对每个参数所有可能组合的一个积分。

它存在的目的仅仅是让整件事加起来等于一。

一个参数时，网格取 100 个点，求值 100 次，咖啡还没落杯就做完了。

五个参数时，是一百亿次求值。

十个参数时——按现代标准只是小模型——你已经来到 10^20，太阳都要对你的时间线有意见了。

这就是维度灾难。它不是能靠工程绕过的性能问题。网格方法在这里死掉。每次都死。

于是你被卡住了：你有一个想从中采样的分布，却连求值都做不到，因为求值需要一个你算不出来的归一化常数。

听起来像是死路一条。

拯救一切的相消

这里的关键操作是：

MCMC 从不问「这一点上的后验概率是多少」。

它只问「这个新点比我当前所在的位置更好还是更差，好多少或差多少」。

这是一个比值。在这个比值中，证据同时出现在分子和分母上。

于是它消去了。

你永远不需要去计算那个不可能算出来的东西，你只需要让自己根本用不着它。

这样做的结果是，你只需要后验分布「差一个常数」的形式，也就是似然乘以先验，而这个你总能算出来。

这就是支撑整个领域的关键洞见。

Metropolis-Hastings 完整算法

最古老的 MCMC 算法出自 Metropolis 等人 1953 年的论文，后来由 Hastings 在 1970 年 推广。它短到可以完全装进脑子里。

你正站在某个参数值上。你想迈出一步。

提出（Propose）。 在当前位置附近抽取一个候选值，通常用以为当前位置中心的高斯分布。

打分（Score）。 计算 R = p(proposed) / p(current)，用未归一化的后验——这是你拥有的全部，也是你需要的全部。

决定（Decide）。 如果 R >= 1，新位置更好。移动过去。总是如此。

如果 R < 1，新位置更差。仍然以概率 R 移动过去。

最后这一行是大家最容易略过的，但正是这一行让算法真正生效。

如果只接受上坡移动，你写出来的就只是一个爬山算法。它会冲刺到最近的峰顶，坐在上面，信心十足地把那个峰当作答案汇报出来，完全没见过分布的其余部分。

偶尔刻意走一步坏棋，链才能从一座山上滚下去，找到另一座山。正是这一步，把贪心优化器变成了采样器。

接受规则不是随意定的。它是被精心构造出来的，目的是让链的平稳分布恰好等于你交给它的后验。

链对高概率区域访问频繁，对低概率区域访问稀少，比例分毫不差。

所以它看起来像随机游走，某种意义上也确实如此——但这是一场被做了手脚的随机游走。

下面是完整算法，货真价实的完整：

import numpy as np

def metropolis(log_post, start, n_steps, step_size):
    x, lp = start, log_post(start)
    chain = []
    for _ in range(n_steps):
        candidate = x + np.random.normal(0, step_size, size=np.shape(x))
        lp_new = log_post(candidate)
        # 对数空间，比值变成减法，什么都不会溢出
        if np.log(np.random.rand()) < lp_new - lp:
            x, lp = candidate, lp_new
        chain.append(x)          # 注意：即使被拒绝也会追加
    return np.array(chain)

里面有两处容易绊倒人的地方。

我们是在对数空间里工作。 在 float64 下，后验很快会下溢成 0，所以比值变成了对数密度的减法。np.log(rand()) < lp_new - lp 与原来的规则完全等价，只是数值上能活下来。

即使被拒绝，也会追加 `x`。 原地不动也是一种真实结果。如果只记录被接受的移动，你会系统性地少算那些尖峰——而尖峰恰恰是大多数提议被拒绝的区域。

画成流程图：

flowchart TD
    A[选一个起始值] --> B[从以它为中心的高斯分布中<br/>提出一个新值]
    B --> C[计算 R = p_new / p_current]
    C --> D{R >= 1?}
    D -- 是，上坡 --> E[接受移动]
    D -- 否，下坡 --> F{掷硬币结果<br/>小于 R?}
    F -- 是 --> E
    F -- 否 --> G[拒绝，原地不动<br/>并再次记录旧值]
    E --> H[记录样本]
    G --> H
    H --> I{仍在燃烧期?}
    I -- 是 --> J[丢弃这个样本]
    I -- 否 --> K[保留它]
    J --> B
    K --> L{样本够了吗?}
    L -- 否 --> B
    L -- 是 --> M[对样本求平均<br/>得到任何你想要的结果]

    classDef decision fill:#f4d35e,stroke:#b8991f,color:#1a1a1a
    classDef start fill:#e9ecef,stroke:#6c757d,color:#1a1a1a
    classDef good fill:#5ee6c8,stroke:#1f9c86,color:#1a1a1a
    classDef bad fill:#ff9a5c,stroke:#c65f22,color:#1a1a1a
    classDef work fill:#6ea8ff,stroke:#2f5fbf,color:#1a1a1a

    class D,F,I,L decision
    class A,M start
    class E,K,H good
    class G,J bad
    class B,C work

步长是唯一的旋钮，它终会反咬你一口

step_size 看起来只是个调参细节。它不是。它决定了链是正常工作，还是对你撒谎。

提议太窄时，几乎所有提议都会被接受，因为你几乎没动。链以冰川般的速度在分布上缓慢爬行，想看任何结果都需要数百万个样本。

提议太宽时，几乎所有提议都会落在糟糕的地方而被拒绝。链动不动就停滞数百步，你的"10,000 个样本"实际上只是大约四十个不同值的重复。

这两种失败从外面看都像是成功。无论哪种情况，你都能拿到一个包含 10,000 个数字的数组。

针对简单随机游走版本的民间经验法则是把接受率控制在 20% 到 25% 左右，这来自 Roberts、Gelman 和 Gilks 一篇确实非常漂亮的渐近工作。打印接受率。永远要打印接受率。

这也是为什么没人会在生产环境手写上面的循环。现代采样器如 Stan 和 PyMC 使用 Hamiltonian Monte Carlo 和 NUTS，它们利用后验的梯度来提出更聪明、更远的移动，而不是盲目的局部抖动，并且会自动调参。核心思想是一样的，只是提议方式不再那么蠢。

燃烧期，或者：扔掉你付过钱的工作

你必须在某个地方启动链，而你的"某个地方"很可能是错的。

如果你的初始猜测落在后验概率很差的区域，链需要花上一段时间走出荒野，才能找到参数空间中真正重要的部分。

这些早期样本是真实样本，消耗了真实的算力，但它们是垃圾。它们描述的是你的糟糕猜测，而不是后验分布。

但请记住马尔可夫性质。链不记得自己从哪里开始。一旦到达高概率区域，它就会留在那里，未来行为不会受到跋涉过程的污染。

所以修复方法简单粗暴得令人尴尬：删掉前几百或前几千个样本。

这就是燃烧期。它不是歪招，而是无记忆性的直接结果。

检查方法就是上图左侧的迹线图：把参数值对迭代次数画出来。健康的链看起来像一条模糊的横向毛毛虫。如果链还在爬升、漂移，或者长时间平直，那它在告诉你一些事情，而且不是什么好消息。

实践中人们会从不同起点跑多条链，检查它们是否收敛到同一处，这正是 R-hat 统计量 所衡量的。

拿到一堆样本之后怎么办

一旦你从后验中拿到了样本，最难的部分就结束了，之后的事情都简单得令人尴尬。

想要参数的均值？对样本求平均。

想要 95% 可信区间？排序后取中间 95%。

想要参数大于 10 的概率？数一数有多少个满足，再除以样本总数。

chain = metropolis(log_post, start=0.0, n_steps=50_000, step_size=0.8)
samples = chain[5_000:]                       # 丢弃燃烧期

samples.mean()                                # 点估计
np.percentile(samples, [2.5, 97.5])           # 95% 可信区间
(samples > 10).mean()                         # 直接得到 P(theta > 10)

你对分布能提出的每个问题都是一个期望，而每个期望都可以用样本均值来近似。这就是蒙特卡洛那一半在最后安静地做着自己的工作。

到这里，名字终于说得通了：马尔可夫链给你样本，蒙特卡洛给你答案，哪一半都无法单独工作。

两半在此汇合

flowchart TD
    P[无法积分的后验] --> Q{能抽取<br/>独立样本吗?}
    Q -- 能 --> MC[普通蒙特卡洛<br/>掷飞镖，取平均]
    Q -- 不能 --> R{至少能求出<br/>差一个常数的值吗?}
    R -- 不能 --> STUCK[彻底卡住了]
    R -- 能 --> CHAIN[构建一条马尔可夫链<br/>其平稳分布<br/>就是那个后验]
    CHAIN --> WALK[让链走很长时间]
    WALK --> S[相关的样本<br/>平均后仍是正确答案]
    MC --> ANS[均值、区间、概率]
    S --> ANS

    classDef decision fill:#f4d35e,stroke:#b8991f,color:#1a1a1a
    classDef start fill:#e9ecef,stroke:#6c757d,color:#1a1a1a
    classDef chip fill:#5ee6c8,stroke:#1f9c86,color:#1a1a1a
    classDef accel fill:#9d8cff,stroke:#5b4bcc,color:#1a1a1a
    classDef bad fill:#ff9a5c,stroke:#c65f22,color:#1a1a1a

    class Q,R decision
    class P start
    class MC,S,ANS chip
    class CHAIN,WALK accel
    class STUCK bad

一个诚实的提醒：我不想把这个方法吹过头。

MCMC 样本是彼此相关的。相邻两步在构造上就注定离得很近，所以 10,000 个 MCMC 样本携带的信息量少于 10,000 个独立样本。真正要紧的指标是有效样本量（effective sample size），它可能比你实际生成的样本数小得多。每个像样的库都会报告这个数字。要去看它。

MCMC 也并非自动生效。它最终会收敛——这句话里「最终」两个字的分量很重。一个具有两个相隔很远的峰、中间夹着深谷的后验分布，能让随机游走链在一个峰上困住比你愿意等待更长的时间，而且整条链全程看上去完全正常。收敛诊断不是被迫害妄想，而是工作本身。

那这跟你有什么关系

因为这个模式完全可以推广到统计学之外。

你有一个无法枚举的对象。你算不出它的总和。但你可以廉价地比较两个候选，也可以随机迈出一步。

这就够了。这就是全部前提。

模拟退火靠的就是这个。统计物理的很大一部分也是。PageRank 的随机冲浪者也是——它是一条马尔可夫链，其平稳分布就是网页排名。当下每一个正在生成图像的扩散模型也是：它们迭代地把噪声引向一个「学到的」分布，而不是你写下来的分布。

1953 年的那篇论文研究的是盒子里的硬球。

后来 SIAM 把 Metropolis 算法列入「二十世纪十大算法」榜单，与快速傅里叶变换（FFT）和 QR 分解并列。

所以下次有人问你：在一个算不出来的分布上怎么做推断？你有答案了。

你不去计算它。

你构建一条链替你在其中行走，让行走替你完成算术。

<em>

你们团队的注意力是有限的，而 AI 生成代码的洪流让「在不拖慢速度的前提下保证生产环境安全可靠」这件事变得越来越难。

我正在构建 LiveReview——一个具备爆炸半径（blast radius）感知能力的 AI 代码评审工具，专为你的关键业务系统打造。

它不会对每个 diff 一视同仁，而是按爆炸半径给每处变更打分——这项变更沿着调用图能波及多远——让你把注意力集中在真正要紧的地方。

把代码评审的精力花在业务风险最高的地方，而不是均匀分摊到每一个 diff 上。

在 GitHub 上为它点个 Star：

{% embed https://github.com/HexmosTech/LiveReview %}

点击下方用你的代码库试用 LiveReview：

</em>

<a href="https://hexmos.com/livereview"><img src="https://dev-to-uploads.s3.us-east-2.amazonaws.com/uploads/articles/vls0pq7nymbrll98je6s.png" alt="LiveReview 横幅" /></a>

原文：https://dev.to/lovestaco/markov-chain-monte-carlo-the-1953-algorithm-hiding-under-modern-ai-5cb4（作者 @lovestaco）