📊 statistics

Hash-augmented adaptive multilevel splitting Monte Carlo algorithm for accurate estimation of two-sample permutation test p-values

本文介绍了一种哈希增强的自适应多层分裂蒙特卡洛算法，该算法已在 Python 包 `hamstest` 中实现，旨在为具有复杂统计量的两样本置换检验精确估计任意小的 p 值，同时解决与分布离散性相关的挑战并确保置信区间有效。

原作者： Nikita Golikov, Vladimir Sukhov, Gennady Korotkevich, Alexey Sergushichev

发布于 2026-07-15

📖 1 分钟阅读☕ 轻松阅读

原作者： Nikita Golikov, Vladimir Sukhov, Gennady Korotkevich, Alexey Sergushichev

原始论文采用 CC BY 4.0 许可（http://creativecommons.org/licenses/by/4.0/）。 ✨ 这是对下方论文的AI生成解释。它不是由作者撰写或认可的。如需技术准确性，请参阅原始论文。阅读完整免责声明

想象一下，你是一名侦探，正试图在一座拥有数百万人口的城市中抓捕一名极其罕见的罪犯。你有一份嫌疑人名单（你的数据），你想知道：“这个特定的线索模式纯粹靠运气发生的可能性有多大？”在统计学世界中，这被称为置换检验（permutation test）。你将线索进行数百万次的随机打乱，以观察这种“幸运”的模式出现的频率。

通常情况下，如果这种模式很常见，你只需要计数那些幸运的打乱情况即可。但如果这种模式极其罕见，比如在万亿次尝试中仅出现一次呢？这就像是在一颗行星大小的沙滩上寻找一颗特定的沙粒。如果你尝试通过随机挑选沙粒的方式（旧有的蒙特卡洛方法）来寻找，你可能会耗费一生去捡沙子，却依然找不到那颗特定的沙粒。为了得到一个像 $10^{-10}$ 这样微小的概率值，你可能需要挑选 $10^{10}$ 颗沙粒，这完全是不切实际的。

问题所在：“卡住”的电梯

这篇论文的作者意识到，当处理这些极小的概率时，标准方法会撞上一堵墙，尤其是因为“沙粒”（数据组合）并不是唯一的。有时，成千上万种不同的打乱方式都会产生完全相同的得分。这就像一部电梯只能停在 1 层、10 层和 100 层，却跳过了 2 到 99 层。如果你想去 99 层，电梯根本无法停在那里，因为它并不存在。这种“离散性”导致数学计算陷入僵局，使得估计一个事件到底有多罕见变得不可能。

解决方案：“标签”与“分层阶梯”

由 Nikita Golikov 领导的团队开发了一个名为 hamstest 的新工具。他们的秘诀在于一个巧妙的技巧，叫做哈希增强自适应多级分裂（hash-augmented adaptive multilevel splitting）。

以下是它的工作原理，使用一个有趣的类比：

阶梯（多级分裂）： 与其尝试直接跳到山顶（稀有事件），不如搭建一个阶梯。他们从底部开始，询问：“有多少人能到达第一级？”然后，“这些人的当中，又有多少人能到达第二级？”他们通过不断将人群拆分为越来越小的群体来向上攀爬。这把一次不可能的飞跃变成了系列易于处理的步骤。
“哈希”标签（修复卡住的电梯）： 大问题在于，许多人都站在同一级台阶上（相同的得分），导致无法进一步拆分人群。为了解决这个问题，作者给每一个人都发了一个独特的、隐形的哈希标签（一个随机数）。即使两个人的得分完全相同，他们的哈希标签也是不同的。这使得算法可以说：“好吧，我们不能按得分来拆分，但我们可以按哈希标签来拆分。”它将一个平坦、卡住的楼层变成了一个平滑、连续的阶梯，使算法总能找到下一个台阶。

他们的发现（以及没发现的）

作者在两个经典的统计检验上测试了这种新方法：Kolmogorov–Smirnov 检验和 Mann–Whitney U 检验。

结果： 在他们的模拟中，这种新方法表现出了惊人的准确性。当他们尝试估计低至 $10^{-243}$ （即 1 后面跟着 243 个零！）的概率时，该方法的估计值精准地落在了真实值上。他们还计算了置信区间（真实答案可能隐藏的范围），在大约 95% 的测试运行中，真实答案都在该范围内。
“全量重采样”规则： 他们尝试了运行模拟的几种不同方式。他们发现，一种叫做**“全量重采样”*的方法（即在每一步都打乱所有*样本）是最可靠且最稳健的。他们建议将 $\alpha = 1$ 这一特定设置作为默认值，因为在他们的测试中效果最好。
他们排除了什么： 他们明确展示了旧的方法（仅使用得分而不使用哈希标签）在数据存在“大幅跳跃”或许多并列值（ties）时会失效。他们证明了如果没有哈希标签，算法会卡住并给出错误的答案。他们还指出，虽然他们的方法对于单侧检验（寻找单一方向的模式）非常有效，但 Kolmogorov–Smirnov 检验的两侧版本比较棘手，因为“电梯”在最顶端可能会发生断层，需要特殊处理。

速度如何？

团队测量了该算法在现代计算机（Apple M3 Pro）上的运行时间。他们发现，所需时间主要取决于事件的稀有程度。如果你正在寻找极其罕见的事件（例如 $p$ 值为 $10^{-100}$ ），它会花费更长时间，因为你必须爬更多的阶梯。然而，对于 Mann–Whitney U 检验，运行时间并不怎么受数据集大小的影响，因为该特定检验的数学计算更新效率非常高。

核心结论

作者并没有“解决”宇宙中所有的统计问题，但他们构建了一个非常强大且灵活的工具，适用于任何科学家可能发明的自定义检验统计量。他们将这个工具封装进了一个名为 hamstest 的免费 Python 库中。

他们建议，对于大多数人来说，使用带有 $\alpha = 1$ 的全量重采样方法是最佳选择。他们还指出，虽然他们的方法很快，但确切的运行时间取决于你所运行的检验的具体数学逻辑。如果你是一名正在处理极小概率和复杂数据的研究人员，这个工具提供了一种无需等待宇宙热寂就能获得准确答案的方法。

简而言之，他们把一部损坏、卡住的电梯变成了一部平滑、高速的自动扶梯，即使路径充满了坑洼，也能带你到达统计学山的顶峰。

技术摘要：用于置换检验 p 值的哈希增强自适应多层分割法

问题陈述
非参数置换检验（如 Kolmogorov–KS 和 Mann–Whitney U 检验）是统计假设检验中的基础方法，特别是在不假设特定数据分布的情况下比较两组数据。虽然精确的 p 值计算是理想的，但由于复杂或自定义检验统计量的可能排列组合数量呈指数级增长（ $\binom{n+m}{n}$ ），这在算法上往往是难以处理的。因此，研究人员依赖蒙特卡洛采样。然而，标准的蒙特卡洛方法在估计极小 p 值（例如 $10^{-10}$ 或更小）时，相对精度较差，为了在多次假设检验校正后达到统计效能，需要极其庞大的样本量。

此外，将多层分割法（一种通过将稀有事件概率分解为条件概率乘积来估计的技术）应用于置换检验时，会引入特定的挑战。置换检验运行在离散空间中，其检验统计量分布通常表现出巨大的跳跃（离散性）。在这种情况下，许多不同的组合可能会产生相同的统计量值。这种离散性可能会阻碍多层分割所需的自适应选择中间层，导致算法失效，或者当无法将移动到下一层的条件概率维持在接近目标常数（例如 0.5）时，产生偏差的估计。

方法论
作者提出了一种哈希增强自适应多层分割蒙特卡洛算法，以解决置换检验中固有的离散性问题。其核心方法论包含三个关键组成部分：

哈希增强排序： 为了解决多个组合共享相同统计量值的问题，作者引入了一种严格的全序关系。他们不仅根据统计量 $S(A)$ 进行排序，还定义了一个基于统计量值首先进行排序，若统计量相等则依据哈希值 $h(A)$ 进行排序的顺序 $A \succ_{S,h} B$ 。
- 哈希函数 $h(A)$ 将组合映射为整数。一种实际的实现方式是使用分配给每个元素的独立随机整数之和，从而允许在马尔可夫链蒙特卡洛（MCMC）步骤中进行 $O(1)$ 的更新。
- 这种增强有效地“平滑”了离散分布，确保对于任何层边界，都能找到一个有效的下一层边界，使得条件概率是可估计的。
结合 MCMC 的自适应多层分割： 该算法通过在由哈希增强顺序定义的嵌套事件序列上的条件概率乘积，来估计 p 值 $Pr(S \ge \gamma)$ 。
- 自适应层： 从一个平凡层开始，算法对组合进行采样，并利用它们的分位数（例如中位数）来定义下一个层边界。
- MCMC 采样： 为了为下一层生成样本，算法采用类 Metropolis 的 MCMC 程序。它从满足新边界条件的现有样本开始，并通过扰动它们（交换元素）来生成新的独立样本。
- 重采样策略： 本文评估了“全重采样”（复制现有的有效样本并扰动所有样本）与“部分重采样”（仅扰动一部分样本）的区别。研究发现，全重采样在应对由问题离散性质引入的采样偏差方面更为稳健。
对数空间估计： 算法通过估计 p 值的对数来处理极端值，以避免数值下溢。它利用 digamma ( $\psi$ ) 和 trigamma ( $\psi_1$ ) 函数来计算对数估计的均值和方差，这些推导源于均匀分布序统计量的性质。这使得构建有效的置信区间成为可能。

主要贡献

算法创新： 引入了哈希增强排序方案，使自适应多层分割能够在离散的置换空间中正常运行，克服了标准方法失效的“大跳跃”问题。
稳健实现： 对重采样策略（全重采样 vs 部分重采样）和 MCMC 停止准则进行了详细分析，证明了采用特定缩放参数（ $\alpha=1$ ）的“全中位数”全重采样在减少偏差和提高收敛速度之间达到了最佳平衡。
软件发布： 开发并开源了 Python 包 hamstest，该包实现了所提出的算法。它允许用户定义自定义检验统计量，并提供了高效的接口（包括 C++ 后端），用于高性能的 p 值估计。
验证： 通过与 Kolmogorov–KS 和 Mann–Whitney U 检验的精确算法进行严格对比，证明了该方法能够准确估计低至 $10^{-243}$ 的 p 值，并提供有效的置信区间。

结果
作者在各种样本量（ $n+m$ 高达 10,000）和目标 p 值（低至 $10^{-100}$ ）下，对 Kolmogorov–KS 和 Mann–Whitney U 检验进行了评估。

准确性： 估计的 p 值与真实值（在有精确方法可用时）高度吻合，没有明显的偏差。在所有测试参数下，估计值的中位数均与真实 p 值一致。
置信区间： 构建的 95% 置信区间（基于 $\log \hat{p} \pm 2\hat{\sigma}$ ）在绝大多数独立运行中成功捕捉到了真实的对数 p 值，验证了理论方差估计的正确性。
性能： 运行时间主要取决于 p 值的大小和样本量。算法表现出相对于所需层数的线性扩展性，使其能够高效估计对于标准蒙特卡洛方法而言在计算上难以实现的极小 p 值。
重采样对比： “全中位数”变体（采用中位数边界选择的全重采样）在收敛速度和对偏差的稳健性方面始终优于部分重采样变体，尤其是在分布离散性更显著的 Kolmogorov–KS 检验中。

意义
本文声称，这项工作是对置换检验中经典蒙特卡洛 p 值估计方法的重大改进。通过解决离散空间中的特定陷阱——即由于统计量值并列而无法定义中间层的问题——所提出的算法能够实现对任意小 p 值的准确且高效的估计。这种能力对于现代高通量应用（如基因组学、神经科学）至关重要，因为在这些领域中，同时进行成千上万次检验，且精确控制错误发现率需要可靠的尾部概率估计。作者指出，该算法的一个先前版本（未进行哈希增强）在基因集富集分析（GSEA）中已经取得了成功；当前的工作将这一效用扩展到了更广泛的两样本统计量类，并解决了在使用特定输入分布时遇到的稳定性问题。hamstest 包的可用性确保了这些先进的统计能力可以被用于实际研究工作流中的用户自定义统计量。