SciPy stats.qmc 准蒙特卡洛（QMC）采样完全指南：从低差异序列到自定义 QMCEngine

本指南以 SciPy 官方教程 doc/source/tutorial/stats/quasi_monte_carlo.rst 为核心骨架，结合 scipy/stats/_qmc.py 源码与 scipy/stats/tests/test_qmc.py 测试用例，系统讲解 scipy.stats.qmc 模块的原理、API 与实战要点。读完本文，你将掌握：为什么 QMC 能对抗"维度灾难"、如何用 discrepancy 量化采样均匀性、如何驱动 Sobol' 与 Halton 序列（含状态续取、重置与快进），以及如何基于 QMCEngine 基类自定义自己的采样引擎。

从 Monte Carlo 到 Quasi-Monte Carlo

Monte Carlo 的基本思想与局限

Monte Carlo（MC）方法是一类依赖重复随机采样来获得数值结果的计算算法，其核心概念是"用随机性去解决原本可能是确定性的问题"。MC 广泛用于物理与数学问题，尤其适合其他方法难以或无法处理的情形，主要覆盖三类问题：最优化、数值积分、从概率分布中抽样。

生成具有特定性质的随机数，比听起来复杂得多。简单 MC 方法生成的样本要求独立同分布（IID），但多次生成随机点集会得到差异巨大的结果。下图（源脚本 doc/source/tutorial/stats/plots/qmc_plot_mc.py）展示了两次完全独立的随机采样——两组点都没有参考先前已生成的点，空间中明显存在未被探索的区域。在仿真中这很危险：一组特定点可能会触发完全不同的行为。

MC 的一大优点是具有已知的收敛性质。以 5 维平方和均值问题为例：

$$f(\mathbf{x}) = \left( \sum_{j=1}^{5}x_j \right)^2,$$

其中 $x_j \sim \mathcal{U}(0,1)$。该函数有解析均值 $\mu = 5/3+5(5-1)/4$。用 MC 数值逼近该均值，近似误差的理论收敛速率为 $O(n^{-1/2})$（对应绘图 doc/source/tutorial/stats/plots/qmc_plot_conv_mc.py）。

虽然收敛有保障，但实际使用者往往希望探索过程更具确定性。普通 MC 可以用随机种子获得可重复的结果，但固定种子会破坏收敛性质：某个种子可能对一类问题有效，对另一类则失效。

维度灾难：网格法的指数爆炸

确定性遍历空间的常见做法是使用覆盖所有参数维度的规则网格（saturated design，饱和设计）。考虑单位超立方体（各维边界为 0 到 1）：若点间距为 0.1，填充单位区间需要 10 个点；2 维超立方体同样间距需要 100 个点；3 维则需要 1,000 个点。随着维度增长，所需实验点数随空间维数呈指数增长，这就是著名的"维度灾难"（curse of dimensionality）。教程给出的示例代码如下（对应绘图 doc/source/tutorial/stats/plots/qmc_plot_curse.py）：

>>> import numpy as np
>>> disc = 10
>>> x1 = np.linspace(0, 1, disc)
>>> x2 = np.linspace(0, 1, disc)
>>> x3 = np.linspace(0, 1, disc)
>>> x1, x2, x3 = np.meshgrid(x1, x2, x3)

QMC：确定性、低差异、可延续

为缓解维度灾难，QMC 方法被设计出来。与 MC 的最大区别在于：QMC 的点不是 IID 的，而是依赖先前生成的点，因此有些方法也被称为"序列"（sequences）。QMC 方法具有三个关键特性：

• 确定性：同一种子/参数下结果可复现；
• 良好的空间覆盖：点分布更均匀，边界附近采样更好，簇与空隙更少；
• 可延续性：部分序列可以继续生成新点而保持良好性质。

下图（doc/source/tutorial/stats/plots/qmc_plot_mc_qmc.py）对比了两组各 256 个点：左侧是普通 MC，右侧是使用 Sobol' 方法的 QMC 设计。可以清楚看到 QMC 版本更均匀，边界采样更好，簇与间隙更少。

评估均匀性的一种方式称为偏差（discrepancy）：Sobol' 点的偏差优于原始 MC。回到均值计算问题，QMC 方法对误差有更好的收敛速率：对该函数可达到 $O(n^{-1})$，对非常光滑的函数甚至更好。下图（doc/source/tutorial/stats/plots/qmc_plot_conv_mc_sobol.py）表明 Sobol' 方法的收敛速率为 $O(n^{-1})$。更多数学细节参见 scipy.stats.qmc 的 API 文档。

用 discrepancy 量化采样均匀性

discrepancy（偏差）是衡量超立方体内一组样本空间填充质量的均匀性准则：它量化了超立方体上的连续均匀分布与 $n$ 个离散样本点构成的经验分布之间的距离。值越低，样本对参数空间的覆盖越好。

考虑两组点，从图上（doc/source/tutorial/stats/plots/qmc_plot_discrepancy.py）明显看出左侧设计比右侧覆盖了更多空间。这可以用 qmc.discrepancy 定量衡量：

>>> import numpy as np
>>> from scipy.stats import qmc
>>> space_1 = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
>>> space_2 = np.array([[1, 5], [2, 4], [3, 3], [4, 2], [5, 1], [6, 6]])
>>> l_bounds = [0.5, 0.5]
>>> u_bounds = [6.5, 6.5]
>>> space_1 = qmc.scale(space_1, l_bounds, u_bounds, reverse=True)
>>> space_2 = qmc.scale(space_2, l_bounds, u_bounds, reverse=True)
>>> qmc.discrepancy(space_1)
0.008142039609053464
>>> qmc.discrepancy(space_2)
0.010456854423869011

注意这里先用 qmc.scale(..., reverse=True) 把位于 <a href="https://link.gitcode.com/i/77aa4f731ed6ac4cb2bdb2c14b30101b" target="_blank">0.5, 6.5] 区间的原始点反向映射回单位超立方体，因为 discrepancy 要求样本落在 $[0,1)^d$ 内——这正是源码中 _ensure_in_unit_hypercube 的硬性校验（[scipy/stats/_qmc.py，超出范围会抛出 ValueError）。

从源码看（scipy/stats/_qmc.py），discrepancy 支持四种方法：

method	全称	含义
CD（默认）	Centered Discrepancy	子空间涉及超立方体一个角
WD	Wrap-around Discrepancy	子空间可环绕边界
MD	Mixture Discrepancy	CD/WD 的混合，覆盖更多准则
L2-star	L2-star Discrepancy	类 CD 但对旋转不变

其中 CD、WD、MD 分别实现文献 2] 的方程 9、10、18（不取平方根），L2-star 计算文献 [3] 方程 10 的量（取平方根）。此外还支持 iterative=True 模式：可计算"好像有 $n+1$ 个样本"时的偏差，配合 update_discrepancy 在逐一评估候选点时比全量重算快得多（对应测试 [scipy/stats/tests/test_qmc.py 中的 test_update_discrepancy）。合理的均匀性度量应满足：因子/行置换不变、坐标旋转不变、可度量低维子投影的均匀性、有几何意义、易计算、满足 Koksma–Hlawka 型不等式、与实验设计其他准则一致。

使用 QMC 引擎：Sobol' 与 Halton

scipy.stats.qmc 实现了多种采样引擎，其中最常用的是 Sobol' 序列与 Halton 序列（对比图见 doc/source/tutorial/stats/plots/qmc_plot_sobol_halton.py，该脚本本身也演示了 qmc.discrepancy 对两类样本的调用）。

必须注意的警告

QMC 方法需要特别小心，使用者必须阅读文档以避免常见陷阱。 例如 Sobol' 要求点数遵循 2 的幂；此外，thinning（抽稀）、burning（烧掉前几个点）或其他点选择操作都会破坏序列的性质，得到的点集可能并不比 MC 好。

这条警告在源码层面得到了印证：Sobol._random 在 num_generated == 0 且 n 不是 2 的幂时会发出警告（scipy/stats/_qmc.py），且序列生成到 $2^{bits}$ 个点后会重复并直接报错；random_base2(m) 则通过 n = 2**m 保证平衡性质，若调用前已生成的点使总数不再是 2 的幂则抛出 ValueError（scipy/stats/_qmc.py）。而 Halton._random 允许从 start_index=self.num_generated 开始任意续取或跳过任意数量的点（scipy/stats/_qmc.py），代价是收敛速率比 Sobol' 慢。

引擎是有状态的：续取、跳过与重置

QMC 引擎是有状态的，可以延续序列、跳过若干点或重置。以 Halton 为例，先取 5 个点，再取第二批 5 个点：

>>> from scipy.stats import qmc
>>> engine = qmc.Halton(d=2)
>>> engine.random(5)
array([[0.22166437, 0.07980522],  # random
       [0.72166437, 0.93165708],
       [0.47166437, 0.41313856],
       [0.97166437, 0.19091633],
       [0.01853937, 0.74647189]])
>>> engine.random(5)
array([[0.51853937, 0.52424967],  # random
       [0.26853937, 0.30202745],
       [0.76853937, 0.857583  ],
       [0.14353937, 0.63536078],
       [0.64353937, 0.01807683]])

重置序列后，再次请求 5 个点会得到与第一次完全相同的 5 个点：

>>> engine.reset()
>>> engine.random(5)
array([[0.22166437, 0.07980522],  # random
       [0.72166437, 0.93165708],
       [0.47166437, 0.41313856],
       [0.97166437, 0.19091633],
       [0.01853937, 0.74647189]])

也可以先快进序列，得到与第二批完全相同的 5 个点：

>>> engine.reset()
>>> engine.fast_forward(5)
>>> engine.random(5)
array([[0.51853937, 0.52424967],  # random
       [0.26853937, 0.30202745],
       [0.76853937, 0.857583  ],
       [0.14353937, 0.63536078],
       [0.64353937, 0.01807683]])

从源码看，random(n, *, workers=1) 每次调用都会累加 self.num_generated（scipy/stats/_qmc.py）；reset() 通过深拷贝保存的 rng_seed 恢复随机数生成器并把 num_generated 归零（scipy/stats/_qmc.py）；fast_forward(n) 则是"生成并丢弃"——内部直接调用 self.random(n=n)（scipy/stats/_qmc.py）。对应测试见 scipy/stats/tests/test_qmc.py 的 test_continuing、test_reset、test_fast_forward。

注意：默认情况下，Sobol 与 Halton 都会进行 scrambling（随机打乱）。打乱后的收敛性质更好，并能防止高维下出现条纹或明显模式。实际上没有理由不使用打乱版本。

从源码确认：Halton 的 scrambling 基于 A. B. Owen, "A randomized Halton algorithm in R", 2017）；Sobol 的 scrambling 采用 LMS+shift（线性矩阵打乱 + 数字随机位移，scipy/stats/_qmc.py），打乱后的 Sobol' 序列适用于奇异被积函数、提供误差估计手段并可提升收敛速率。

两大引擎的关键参数

Sobol'（scipy/stats/_qmc.py）：

• d：序列维度，最大维度为 21201（MAXDIM，方向数来自 Joe & Kuo 2008 的预计算表）；
• scramble：默认 True，使用 LMS+shift 打乱；
• bits：控制可生成的最大点数 $2^{bits}$，默认 None（兼容旧版本时为 30），最大 64；超过上限会抛 ValueError；
• 关键方法 random_base2(m)：安全地抽取 $n=2^m$ 个点，保证序列的平衡性质；
• optimization：可选 None / "random-cd" / "lloyd"（版本 1.10.0 起）。

Halton（scipy/stats/_qmc.py）：

• d：维度，第 $n$ 维使用第 $n$ 个素数作为基的 Van der Corput 序列；
• scramble：默认 True；源码注释指出即使中等维度，Halton 序列也有严重的条纹伪影，打乱可显著改善，并支持基于复制的误差估计与无界被积函数；
• 支持任意数量的点采样或跳过（random / fast_forward），并支持 workers 并行（n > 10^3 时多线程才更快，见 scipy/stats/_qmc.py）。

补充说明：QMCEngine 还提供了 integers(l_bounds, u_bounds, n, endpoint) 方法，可按 $\lfloor (b-a)\cdot sample + a \rfloor$ 将 $<a href="https://link.gitcode.com/i/5fe27355821724611ef3716d1ee7ea5d" target="_blank">0,1)$ 样本映射为任意整数区间，且保持无偏性与方差缩减性质（[scipy/stats/_qmc.py）。

自定义 QMCEngine：子类化指南

要创建自己的 QMCEngine，需要定义几个方法。官方教程给出的示例是包装 numpy.random.Generator：

>>> import numpy as np
>>> from scipy.stats import qmc
>>> class RandomEngine(qmc.QMCEngine):
...     def __init__(self, d, seed=None):
...         super().__init__(d=d, seed=seed)
...         self.rng = np.random.default_rng(self.rng_seed)
...
...
...     def _random(self, n=1, *, workers=1):
...         return self.rng.random((n, self.d))
...
...
...     def reset(self):
...         self.rng = np.random.default_rng(self.rng_seed)
...         self.num_generated = 0
...         return self
...
...
...     def fast_forward(self, n):
...         self.random(n)
...         return self

然后像使用任何其他 QMC 引擎一样使用它：

>>> engine = RandomEngine(2)
>>> engine.random(5)
array([[0.22733602, 0.31675834],  # random
       [0.79736546, 0.67625467],
       [0.39110955, 0.33281393],
       [0.59830875, 0.18673419],
       [0.67275604, 0.94180287]])
>>> engine.reset()
>>> engine.random(5)
array([[0.22733602, 0.31675834],  # random
       [0.79736546, 0.67625467],
       [0.39110955, 0.33281393],
       [0.59830875, 0.18673419],
       [0.67275604, 0.94180287]])

从基类源码（scipy/stats/_qmc.py）可以提炼出子类化契约：

• 必须定义：__init__(d, rng=None)（至少固定维度，确定性方法如 Halton 可省略 rng 参数）、_random(n, *, workers=1)（从引擎抽取 $n$ 个点，workers 用于并行，参考 Halton 的实现）；
• 可选覆写：reset()（将引擎重置到初始状态）、fast_forward(n)（确定性序列如 Halton 可跳过前 $n$ 次抽取）；
• 基类 _initialize 会统一校验维度（d 必须为非负整数）、实例化并深拷贝 rng_seed、初始化 num_generated 计数器，并处理 optimization 配置（random-cd 的 n_nochange=100、n_iters=10000，lloyd 的 tol=1e-5、maxiter=10 等默认参数）；
• 约定：样本分布在半开区间 $[0, 1)$ 内，实例可访问 d 与 rng 属性；
• 需要注意当前仓库处于 SPEC-007 迁移期：构造参数已从 seed 过渡到 rng，过渡期内两者均可使用但只能指定其一，seed 会逐步进入弃用流程。

测试文件 scipy/stats/tests/test_qmc.py 的 test_subclassing_QMCEngine 也验证了子类化后 random、reset、fast_forward 的行为符合预期。

QMC 使用准则（官方要点）

教程在结尾给出了一组必须遵守的使用准则，直接决定你能否获得优于 MC 的收益：

1. QMC 有规则！ 务必阅读文档，否则可能得不到任何优于 MC 的好处；
2. 需要恰好 $2^m$ 个点时使用 Sobol；
3. Halton 可以采样或跳过任意数量的点，代价是收敛速率比 Sobol' 慢；
4. 永远不要移除序列开头的点——这会毁掉序列的性质；
5. Scrambling（打乱）总是更好；
6. 如果使用基于 LHS 的方法，无法在不失去 LHS 性质的前提下增加点数（存在一些扩展方法，但当前未实现）。

实践建议与常见陷阱

结合教程与源码，给出几项落地建议：

• 选型速查：点数必须是 2 的幂、追求 $O(n^{-1})$ 收敛 → Sobol + random_base2；需要灵活点数或任意跳取 → Halton；需要分层但允许随机化 → LatinHypercube（默认 scramble=True，支持 strength=2 的正交阵列 LHS）；更多场景可关注 PoissonDisk、MultivariateNormalQMC、MultinomialQMC 等引擎（scipy/stats/_qmc.py）。
• 缩放与反缩放：所有引擎默认输出 $<a href="https://link.gitcode.com/i/a06ccb71468f21d1b95b305eb49fb8b2" target="_blank">0,1)^d$，用 qmc.scale(sample, l_bounds, u_bounds) 线性映射到真实边界，reverse=True 可做逆变换（[scipy/stats/_qmc.py）。
• 后处理优化：1.10.0 起可在构造时指定 optimization="random-cd"（通过随机坐标置换降低 centered discrepancy）或 "lloyd"（Lloyd-Max 扰动使样本趋于等间距），但需注意这是后处理步骤，不保证保留序列全部性质。
• 质量评估：始终用 qmc.discrepancy（或高维下的 geometric_discrepancy）量化样本质量；评估新候选点时用 iterative=True + update_discrepancy 避免全量重算。

参考资料与仓库位置

• 教程原文：doc/source/tutorial/stats/quasi_monte_carlo.rst
• 核心实现：scipy/stats/_qmc.py（含 scale、discrepancy、QMCEngine、Halton、Sobol 等全部引擎）
• 公开 API 入口：scipy/stats/qmc.py
• 测试用例：scipy/stats/tests/test_qmc.py
• 教程配图脚本：doc/source/tutorial/stats/plots/（qmc_plot_mc.py、qmc_plot_conv_mc.py、qmc_plot_curse.py、qmc_plot_mc_qmc.py、qmc_plot_conv_mc_sobol.py、qmc_plot_discrepancy.py、qmc_plot_sobol_halton.py）