接受-拒绝采样法

玻尔百科

核心要点

接受-拒绝采样法通过从一个更简单的分布中提出候选样本，并根据概率测试接受它们，从而从一个复杂的目标分布中生成样本。
算法的效率与包络常数 $M$ 成反比，通过选择一个与目标分布紧密相似的提议分布可以最小化 $M$ 。
在高维空间中，由于测度集中现象（即所谓的维度灾难），该方法的效率可能会呈指数级下降。
应用范围从简单的几何采样到高能物理中的复杂模拟，以及通过稀疏法对非齐次泊松过程进行建模。

引言

生成遵循特定、通常是复杂概率分布的随机数，是贯穿科学与工程的一项基本任务。虽然计算机可以轻松生成均匀分布的随机数，但挑战在于如何将这种简单的输出转换为能够模拟从金融市场波动到亚原子粒子行为等复杂模式的样本。接受-拒绝采样法为这个问题提供了一个 brilliantly intuitive 和精确的解决方案，是蒙特卡洛模拟的基石之一。本文将深入探讨这一强大的技术。首先，在 原理与机制 部分，我们将通过一个简单的飞镖靶类比来探索该方法的核心逻辑，将其算法形式化，并分析决定其效率的关键因素，包括其在高维空间中的局限性。随后，应用与跨学科联系 一章将展示该方法非凡的通用性，介绍其在从简单几何学到现代高能物理核心的大规模模拟等领域中的应用。

原理与机制

我们如何能让计算机生成遵循非常特定、甚至可能很复杂的模式的随机数？假设我们有一台机器，它只能吐出完全随机、均匀分布的数字——就像掷一个完美平衡、有无限个面的骰子。我们的任务是使用这个简单的工具来生成符合某个更复杂的概率分布的数字，比如说宇宙射线的能量分布或金融市场的波动。这是科学和工程中的一个核心问题，而最优雅的解决方案之一是一个名为 接受-拒绝采样法 的绝妙直观的想法。

飛鏢靶類比：將採樣視為面積問題

想象一下你想要采样的概率分布，即你的 目标分布 $f(x)$ ，就像画在一张纸上的一条曲线。这条曲线下的总面积恰好为1，代表100%的概率。我们的目标是以某种方式生成随机点 $x$ ，使得沿水平轴的点密度与曲线 $f(x)$ 的高度相匹配。

如果我们技术高超，可以发明一种特殊的飞镖投掷机，能以恰到好处的模式将飞镖投掷到这条曲线下方。但这听起来很复杂。让我们试试一种更简单的方法。如果我们画一个能完全包住我们目标曲线的简单矩形呢？我们已经知道如何在矩形内均匀地生成点——这很简单！我们只需随机选择一个水平位置和一个垂直位置。

这就是接受-拒绝采样法的精髓。我们使用一个易于采样的 提议分布，我们称之为 $g(x)$ ，作为我们“矩形”的底部。为确保矩形足够高，能在任何地方都包含我们的目标曲线 $f(x)$ ，我们将 $g(x)$ 乘以一个常数 $M$ 。这个常数，即 包络常数，必须选择得当，使得不等式 $f(x) \le M g(x)$ 对所有可能的 $x$ 值都成立。我们的“矩形”现在是由曲线 $M g(x)$ 定义的形状。

现在，我们来玩一个飞镖游戏。我们投掷的飞镖均匀分布在由 $M g(x)$ 定义的更大形状内。如果飞镖落在我们的目标曲线 $f(x)$ 下方，我们就“接受”它的水平位置作为一个有效样本。如果它落在 $f(x)$ 上方但在矩形内部，我们就“拒绝”它并丢弃。通过重复这个过程，收集到的被接受的水平位置将完全按照我们原始的目标分布 $f(x)$ 分布。

算法实践：如何说“是”或“否”

让我们把这个视觉类比转换成一个精确的算法。这个过程有两个简单的步骤：提议和决策。

提议： 我们从简单的提议分布 $g(x)$ 中抽取一个候选值，称之为 $Y$ 。这就像选择我们飞镖的水平坐标。例如，如果我们想在区间 $[0, 1]$ 上从一个复杂的曲线进行采样，我们可能会选择最简单的提议分布：在该区间上的均匀分布，其中 $g(x)=1$ 。
决策： 为了决定是接受还是拒绝 $Y$ ，我们需要模拟飞镖的垂直位置。我们生成第二个随机数 $U$ ，这次是从0到1的均匀分布中抽取。这个 $U$ 代表飞镖在位置 $Y$ 处“矩形”内的相对高度。在 $Y$ 处的矩形总高度是 $M g(Y)$ 。如果飞鏢的隨機高度（由 $U$ 表示）位於對應於目標的較低部分，則它被「接受」。在数学上，如果我们的随机分数 $U$ 乘以总高度小于或等于目标的高度，我们就接受 $Y$ ：
$U \times (\text{在 } Y \text{ 处的总高度}) \le \text{在 } Y \text{ 处的目标高度}$ $U \cdot M g(Y) \le f(Y)$
稍作整理，我们便得到了著名的 接受条件：
$U \le \frac{f(Y)}{M g(Y)}$

但是，这个异常简单的程序为什么能奏效呢？这感觉像是魔法，但其实是纯粹的逻辑。提议某个特定值 $Y=y$ 的概率由 $g(y)$ 给出。然后接受它的概率是我们的均匀随机数 $U$ 小于比率 $\frac{f(y)}{M g(y)}$ 的概率。由于 $U$ 是均匀的，这个概率就是该比率本身。因此，提议 $y$ 并且它被接受的组合概率是：

P(\text{提议 } y \text{ 并接受}) = g(y) \times \frac{f(y)}{M g(y)} = \frac{f(y)}{M}

这是一个优美的结果！被接受样本的概率密度与我们的目标函数 $f(y)$ 成正比。常数 $M$ 和提议函数 $g(y)$ 从分布的最终形式中消失了，留下了一个纯粹、未失真的目标分布版本。当我们审视所有被接受样本的集合时，它们的分布不仅仅是一个近似——它完全是目标分布 $f(x)$ 。通过这种方法生成的每个样本，如果基于真正独立的随机数，都是从目标分布中进行的一次独立且完美的抽取。

简单的代价：效率与提议的艺术

这个优雅的方法并非没有代价：被拒绝的样本。实际上，我们浪费了一部分计算精力。算法的效率就是任何给定提议被接受的概率。我们可以通过对所有可能的提议的接受概率进行求和（或积分）来计算这个值。

P(\text{接受}) = \int \frac{f(x)}{M} dx = \frac{1}{M} \int f(x) dx

由于 $f(x)$ 是一个概率密度，它在整个空间上的积分是1。因此，总体的 接受概率 就是 $1/M$ 。因此，为了获得一个被接受的样本，我们平均需要生成的提议数量是 $M$ 。

这就揭示了该方法的核心挑战：为了使算法高效，我们必须选择我们的提议分布 $g(x)$ 和包络常数 $M$ ，使得 $M$ 尽可能小。 $M$ 的最小可能值是比率 $\frac{f(x)}{g(x)}$ 的最大值（或上确界）。

M = \sup_{x} \frac{f(x)}{g(x)}

因此，选择一个好的提议分布 $g(x)$ 是一门艺术。它应该易于采样，但其形状应尽可能地模仿目标 $f(x)$ 。如果 $g(x)$ 的形状与 $f(x)$ 非常不同，那么在某些区域比率会非常大，迫使 $M$ 变大，从而使算法效率低下。我们甚至可以更进一步，数学上优化一整族提议分布的参数，以找到最小化 $M$ 从而最大化效率的那个。

此外，高 $M$ 值的代价不仅仅是平均运行时间变慢。每次接受所需的试验次数遵循几何分布。该分布的方差是 $M^2 - M$ 。对于一个大的 $M$ ，运行时间的标准差近似为 $M$ 本身。这意味着，如果你的算法平均效率为1%（ $M=100$ ），那么每个样本的运行时间将会剧烈波动，使得你的模拟性能高度不可预测——这对于任何实际应用来说都是一场噩梦。

维度灾难：当直觉失效时

飞镖靶的类比在一维或二维空间中非常强大。但当我们需要在有成百上千个维度的空间中采样一个点时（这在统计力学或粒子物理学中很常见），情况又如何呢？在这里，我们的低维直觉会 spectacularly地崩溃。

考虑在 $d$ 维空间中从一个标准的多维钟形曲线（高斯分布）进行采样。让我们使用一个稍宽的钟形曲线作为我们的提议，这似乎是一个合理的选择。在一维中，这行得通。但随着维度 $d$ 的增加，惊人的事情发生了。接受概率急剧下降，其缩放关系为 $\sigma^{-d}$ ，其中 $\sigma>1$ 是衡量我们提议分布宽度的指标。效率不仅仅是降低，而是呈指数级地趋向于零。

这是 维度灾难 的一个典型表现。其原因是一个关于高维空间的奇怪几何事实，称为 测度集中现象。在高维空间中，一个球体的几乎所有体积都不在中心附近，而是集中在靠近表面的一个非常薄的“壳”中。我们的目标分布的概率集中在某个特定半径的壳层上，而我们稍宽的提议分布的概率则集中在另一个半径更大的不同壳层上。这两个壳层几乎不重叠。我们几乎总是在目标概率实际上为零的区域提出候选样本。我们的飞镖靶实际上已经变成了全是“拒绝”区域，而“接受”区域则缩小到一个无穷小的斑点。

摆脱这个陷阱的方法是放弃简单、单一的提议分布，转而利用问题的结构。如果目标分布的维度是独立的，我们可以构建一个同样是独立分布乘积的提议分布，这样可以将效率恢复到100%。这一洞见为更先进的蒙特卡洛方法铺平了道路，这些方法旨在驾驭高维空间的险恶地貌。

一点警示：机器中的幽灵

接受-拒绝采样法的整个数学大厦，及其精确性的证明，都建立在一个关键的假设之上：我们有一个完美的均匀随机数来源，特别是对于关键的决策步骤 $U \le \frac{f(Y)}{M g(Y)}$ 。

如果我们的“均匀”随机数生成器有缺陷会怎样？例如，如果它有轻微的偏差，产生较小数的频率比较大数高怎么办？这个看似微不足道的瑕疵可能会产生灾难性的后果。如果 $U$ 的值系统性地偏小，我们的接受条件会比应有的更频繁地得到满足，特别是在阈值 $\frac{f(Y)}{M g(Y)}$ 较小的区域。这将系统性地扭曲被接受样本的分布。最终得到的数据将不再遵循目标分布 $f(x)$ ，而是一个被篡改了的版本。那个优美的正确性数学保证也就烟消云散了。

这给我们上了深刻的最后一课。一个算法的理论优雅性，仅与其实际（或计算）实现的好坏相当。接受-拒绝采样法为到达期望的目的地提供了一张完美的地图，但我们必须相信我们的交通工具——随机数生成器——能够忠实地遵循它。这个基础工具中的任何缺陷，都像机器中的幽灵一样，默默地将我们的模拟引向歧途。

应用与跨学科联系

在理解了接受-拒绝采样法的原理之后，我们可能会倾向于认为它只是一个聪明但或许小众的数学技巧。事实远非如此。这个“提议与检验”的简单思想是计算科学家工具库中最通用、最强大的工具之一。它的应用不局限于单一领域，而是横跨一个广阔的领域，从绘制简单图形到模拟宇宙的构造。让我们踏上一段旅程，看看这一个优雅的原则如何为千差万别的问题提供统一的解决方案。

几何学家的飞镖靶

从本质上讲，接受-拒绝采样法是一个用完美的准头，但在形状奇特的靶子上玩飞镖的游戏。想象一下，你想在一个奇特的形状内均匀地生成点，比如说由抛物线 $y=x^2$ 和直线 $y=1$ 所包围的区域。你会怎么做？一个非常简单的方法是将这个形状包围在一个简单的矩形中——一个我们确实知道如何采样的形状，就像向一个矩形靶子投掷飞镖一样。然后我们向矩形均匀地投掷飞镖。如果飞镖落在我们的目标抛物线区域内，我们就保留它。如果落在外面，我们就丢弃它。就是这样！我们保留的点集将在所期望的抛物线区域内完美、均匀地分布。

这种简单的代价是什么？是我们必须扔掉的飞镖数量。接受任何给定飞镖的概率就是目标区域与提议区域的面积之比。在我们的抛物线例子中，它占据了包围矩形的三分之二，我们预计平均每三个提议中会接受两个。这种效率，即面积之比，是一个非常直观的概念，构成了该方法的基础。

但世界很少是均匀的。我们常常需要从某些结果比其他结果更可能的分布中采样。想象我们的飞镖靶不再是一个简单的形状，而是一个有山丘和山谷的地形，在点 $(x,y)$ 处着陆的概率由某个函数 $f(x,y)$ 给出。我们可以调整我们的游戏。我们仍然使用一个简单的提议区域，比如一个正方形，但在我们的测试中增加第二步。在提议一个候选点 $(X_c, Y_c)$ 后，我们不只是检查它是否在正确的区域内；我们玩一个机会游戏，其成功概率与该点的目标密度 $f(X_c, Y_c)$ 成正比。我们可以将其想象为画出这个点，然后掷一个骰子来决定是否保留它，其中该点地形的“高度”会影响骰子的结果。这样，我们更有可能接受来自“山丘”（高概率区域）的点，而不太可能接受来自“山谷”（低概率区域）的点。这个扩展将我们从简单的几何采样带到了几乎可以想象到的任何概率分布的采样。

超越几何：滴答的时钟和抛掷的硬币

接受-拒绝采样法的威力在于它不仅限于连续的几何空间。它对离散事件同样有效。假设我们想从一个奇特的集合，比如 {1, 2, 3, 4}，模拟一个随机整数，其概率不均匀，而是遵循像 $(1/2)^k$ 这样的模式。我们可以从集合中均匀地提议一个候选者——就像掷一个公正的四面骰子——然后使用我们的接受测试，根据目标概率加权，来决定是否保留结果。基本逻辑保持不变，展示了该方法的深刻普适性。

这种普适性延伸到随时间展开的过程中。考虑对网络网关处数据包到达的情况进行建模。到达率可能不是恒定的；它可能在高峰时段激增，在夜间下降。这是一个非齐次泊松过程。我们如何模拟它？其中一种最优雅的技术，称为“稀疏法”，其实就是接受-拒绝采样法的变相。我们首先从一个更简单的、齐次泊松过程中生成候选事件，其速率 $\lambda_{\max}$ 保证至少与任何时间的真实速率 $\lambda(t)$ 一样高。这是我们的提议分布。然后，对于每个在时间 $t_p$ 提议的事件，我们通过以概率 $\lambda(t_p) / \lambda_{\max}$ 接受它来“稀疏”这个事件流。在高实际流量时段的事件很可能被保留，而在低谷时段的事件很可能被丢弃。最终得到的被接受事件流完美地模拟了我们想要建模的复杂、时变过程。从网络流量到放射性衰变，这一原理为模拟提供了一条直接而直观的途径。

虚拟宇宙：接受-拒绝采样法在高能物理中的应用

接受-拒绝采样法最引人注目的应用可能是在高能物理（HEP）领域，它是解释来自像LHC这样的粒子对撞机数据的海量模拟的基石。当粒子碰撞时，它们可以产生一连串的新粒子，量子力学定律决定了任何特定结果的概率。这个概率由一个看起来很吓人的对象——矩阵元的平方 $|\mathcal{M}|^2$ ——所支配，它是末态粒子动量的函数。

为了检验一个新理论，物理学家必须将其预测与实验数据进行比较。这意味着他们需要生成数百万或数十亿个模拟“事件”，每个事件都是一组根据 $|\mathcal{M}|^2$ 指定的分布抽取的粒子动量。这是一个巨大的采样问题。一个典型的策略是使用一个巧妙的算法来生成满足守恒定律但其他方面是“通用”的动量配置——这是提议分布，由所谓的雅可比因子 $J$ 控制。然后接受-拒绝采样法登场。对于每个提议的事件，计算一个权重 $w = |\mathcalM}|^2 / J$ 。这个权重是真实物理概率与提议概率之比。通过以与此权重成比例的概率接受事件，物理学家产生了一个“无权重”的事件样本，其行为与理论预测完全一致。他们实际上在计算机中创造了一个虚拟宇宙。

在处理真实实验的实际问题时，该方法的优雅之处更加闪耀。物理探测器无法看到所有方向或所有能量的粒子。数据总是有“截断”，例如，要求粒子的横向动量 $p_T$ 高于某个阈值。在模拟中如何处理这些截断？接受-拒绝框架免费处理了它们。对于任何会落在截断之外的事件，目标密度被简单地定义为零。当进行接受测试时，这些事件会自动且确定性地被拒绝 [@problemid:3512604]。采样算法自然地只关注那些对实验具有物理意义的事件。

值得注意的是，这些模拟的效率可以直接与 underlying physics 联系起来。在CP破坏——物质与反物质之间的一种微妙不对稱性——的研究中，衰变产物的分布可能具有像 $1 + \alpha \cos(2\phi + \delta)$ 这样的调制。参数 $\alpha$ 衡量了这种破坏的强度。当使用简单的均匀提议来采样这个分布时，平均接受概率结果为 $p = 1 / (1 + |\alpha|)$ 。这个优美的结果意味着，物理现象越复杂（CP破坏 $|\alpha|$ 越大），模拟的效率就越低，研究它所需的计算精力就越多。发现的成本被写进了算法本身。

算法的艺术：效率、组合与双尾传说

接受-拒绝采样法不仅仅是一个独立的算法；它是一个基本的构建模块。假设需要从一个非常复杂的“混合”分布中采样，该分布既有连续部分也有离散尖峰（原子）。一个强大的策略是分解问题。可以构建一个两阶段采样器，首先决定是从连续部分生成样本还是从离散原子中生成样本，然后对所选组件使用量身定制的接受-拒绝步骤。这展示了组合简单思想所带来的模块化和强大功能。

然而，一位经验丰富的计算科学家知道，拥有一个工具是不够的；必须知道何时使用它。如果我们的目标分布是一个混合体， $p(x) = \sum \pi_k f_k(x)$ ，并且我们可以轻松地从每个分量 $f_k$ 中采样，我们可以使用直接复合方法。或者，我们可以对总密度 $p(x)$ 使用接受-拒绝采样法。哪个更好？答案在于成本效益分析。通过对每个操作的计算成本——从提议中抽取、评估密度等——进行建模，可以得出一个精确的阈值 $M^{\star}$ 。如果接受-拒绝采样法的拒绝常数 $M$ 大于此阈值，那么复合方法预计会更经济。这是工程思维的精髓：根据对权衡的量化理解来选择最适合工作的工具。

这给我们带来了最后一个关键的教训：接受-拒绝采样法的好坏取决于其提议分布。它的阿喀琉斯之踵是要求包络常数 $M = \sup_x p(x)/q(x)$ 必须是有限的。考虑尝试使用像高斯分布这样的轻尾提议从一个重尾分布（衰减缓慢的分布，如幂律分布）中采样。比率 $p(x)/q(x)$ 将在尾部爆炸，因为高斯提议密度比目标密度消失得快得多。常数 $M$ 将是无限的，方法完全失败。这就像试图用一个满是小孔的网来捕鲸。在这种情况下，其他方法如自归一化重要性采样（SNIS）可能仍然有效，即使它们也带来了自己的一系列挑战，比如高方差。

这个警示性的故事揭示了蒙特卡洛模拟的真正艺术。成功往往取决于选择一个“尊重”目标 $p(x)$ 的提议分布 $q(x)$ ，特别是其尾部行为。我们开始时玩的那个简单的飞镖游戏，最终让我们深刻体会到概率分布之间微妙的相互作用。接受-拒绝采样法，无论在其成功还是潜在的失败中，都教给我们计算科学的一个基本原则：要有效地模拟一个现象，你必须首先理解它的本质特征。