SciPy 准蒙特卡洛（scipy.stats.qmc）实战指南：从低差异序列到随机化 QMC 积分

本篇指南围绕 SciPy 官方参考文档中的 scipy.stats.qmc 模块 展开，系统讲解准蒙特卡洛（Quasi-Monte Carlo, QMC）方法的数学背景、scipy.stats.qmc 提供的低差异采样引擎与辅助函数，并结合仓库源码 scipy/stats/_qmc.py 与测试 scipy/stats/tests/test_qmc.py 给出可运行、可验证的实战示例。读完本文你将掌握：如何使用 Sobol、Halton、LatinHypercube、PoissonDisk 等引擎生成低差异点集，如何用 discrepancy / geometric_discrepancy 量化点集均匀性，以及如何将 QMC 点集缩放、旋转、用于多元正态与多项分布采样，从而把数值积分误差从普通 MC 的 O(n^{-1/2}) 推向接近 O(n^{-1})。

一、什么是准蒙特卡洛：模块核心概念

scipy.stats.qmc 模块提供准蒙特卡洛生成器及配套辅助函数。模块的完整介绍见 scipy/stats/qmc.py 的模块级 docstring（由参考页 doc/source/reference/stats.qmc.rst 通过 automodule 指令渲染）。

QMC 方法生成一个 n × d 的、取值于 [0,1] 的超立方体内的数组，可以替代从 U[0,1]^d 分布独立抽取的 n 个随机点。与随机点相比，QMC 点被设计为更少出现空隙（gaps）和聚集（clumps），这一性质由 discrepancy（差异度/偏差） 度量来量化。根据 Koksma-Hlawka 不等式，低差异能降低积分误差的上界：对表现良好的被积函数 f 在 n 个 QMC 点上取平均，积分误差可接近 O(n^{-1})，远优于蒙特卡洛方法的 O(n^{-1/2})。

理解 QMC 需要把握几个关键特性：

• 特殊样本量（special sample sizes）：大多数 QMC 构造是为特定的 n 设计的，例如 2 的幂或较大的素数。即使只把样本量改变 1，也可能显著降低其性能甚至收敛速率。例如，如果某方法是为 n = 2^m 设计的，那么 n = 100 个点的精度可能反而不如 n = 64。
• 可扩展性（extensibility）：有些构造在 n 方向上可扩展——总能找到更大的特殊样本量 n' > n，且往往存在无限递增的特殊样本量序列；有些构造在 d 方向上可扩展，通常不需要 d 取特殊值；还有些构造在 n 和 d 两个方向都可扩展。
• 确定性（deterministic）：QMC 点是确定性的，这导致难以估计积分精度，由此引出随机化 QMC（RQMC）。

1.1 随机化 QMC（RQMC）与 scrambled nets

RQMC 点被构造为：每个点单独来看服从 U[0,1]^d，而整体上 n 个点仍保持低差异。可以通过生成 R 次独立的 RQMC 重复来观察计算的稳定性，再从 R 个独立取值出发，用 t 检验（或 bootstrap t 检验）给出均值的近似置信区间。部分 RQMC 方法能达到 o(1/n) 的均方根误差（RMSE），优于未随机化的 QMC——直观解释是误差由许多小误差求和，随机误差会相互抵消，而确定性误差不会。RQMC 对被积函数奇异、或出于其他原因不满足黎曼可积的情形也更有优势。

Scrambled nets（打乱网） 是一类具备宝贵稳健性质的 RQMC：

• 若被积函数平方可积，其方差满足 var_SNET = o(1/n)；
• 存在一个对所有平方可积被积函数同时成立的有限上界 var_SNET / var_MC；
• 对 p > 1 的 L^p 空间中的 f，满足强大数定律；
• 在若干特殊情形下存在中心极限定理；
• 对足够光滑的被积函数，RMSE 可接近 O(n^{-3})。

1.2 主要构造：lattice rules 与 digital nets

QMC 方法的两大类是晶格规则（lattice rules）和数字网/序列（digital nets and sequences），二者在**多项式晶格规则（polynomial lattice rules）**中交汇——后者可以产生数字网。晶格规则需要搜索好的构造参数，而数字网有被广泛使用的默认构造。

模块文档对主流方法给出的定位如下：

方法	类型	特性
Sobol' 序列	数字网	最广泛使用；在 n、d 两个方向均可扩展；可打乱；特殊样本量为 2 的幂
Halton 序列	类数字网	前几个维度均分性质远好于后面的维度；几乎无特殊样本量；精度通常不及 Sobol'；可打乱
Faure 网	数字网	所有维度等优，但特殊样本量随维度 d 快速增大；可打乱
Niederreiter-Xing 网	数字网	渐近性质最优，但实证表现并不突出
高阶数字网（higher order digital nets）	数字网	通过数字交错过程构造，在被积函数更光滑时可达到更高渐近精度，可打乱；目前缺少实证工作

高阶数字网通过在被构造点的数字中做交错（digit interleaving）来形成，可在更高的光滑性条件下获得更高的渐近精度，并可打乱。

使用 QMC 就像使用一个小型随机数生成器的整个周期——构造类似，因此计算成本也类似。此外：

• Baker 变换（tent 函数）：形式为 1 - 2|x - 1/2|，x 从 0 到 1 时该函数先由 0 升到 1 再降回 0。它非常适合为晶格规则构造周期函数，有时还能改善收敛速率。
• QMC 与 MCMC：把 MCMC 看作在 [0,1]^d 中使用 n = 1 个点、d 极大，可借助相关论文中的提议方法在强条件下获得改进的收敛速率。
• 维度诅咒：(R)QMC 无法战胜 Bahkvalov 维度诅咒——对任意随机或确定方法，都存在使它在高维下表现很差的最坏情形函数（例如在所有 n 个点上为 0、其余地方很大的函数）。当被积函数不是最坏情形时，(R)QMC 相对 MC 能带来巨大改进，尤其适合能由“少量输入变量的函数之和”良好逼近的被积函数。
• 光滑性要求：要看到相对 IID MC 的改进，(R)QMC 需要被积函数具有一定的光滑性——混合一阶偏导 ∂^d f/∂x_1⋯∂x_d 大致需要可积。例如“单位超球内为 1、球外为 0”的指示函数在任意 d = 2 时都有无穷大的 Hardy-Krause 变差。

1.3 Sobol' 的方向数与维度上限

Sobol' 点有很多版本，取决于所谓的方向数（direction numbers）——它们是搜索得到的并已被制表。模块文档特别指出，最广为使用的一组方向数来自 Joe 与 Kuo 的工作，该组方向数支持在维度上扩展到 d = 21201。这一事实在源码 scipy/stats/_qmc.py 的 Sobol 类实现中亦有对应（_qmc_cy 与方向数生成相关代码）。

二、模块 API 全景

scipy.stats.qmc 的公开 API 由 scipy/stats/_qmc.py 中的 __all__ 声明（第 44-46 行）：

__all__ = ['scale', 'discrepancy', 'geometric_discrepancy', 'update_discrepancy',
           'QMCEngine', 'Sobol', 'Halton', 'LatinHypercube', 'PoissonDisk',
           'MultinomialQMC', 'MultivariateNormalQMC']

按参考文档的组织方式可分为两组：

Engines（引擎）

类	说明
QMCEngine	所有 QMC 引擎的抽象基类（ABC），定义统一的采样接口
Sobol	Sobol' 低差异序列，n、d 双向可扩展，支持打乱
Halton	Halton 低差异序列，按不同素数基构造
LatinHypercube	拉丁超立方采样（LHS），支持随机化与优化
PoissonDisk	泊松盘采样，保证点间最小距离
MultinomialQMC	基于 QMC 的多项分布随机数
MultivariateNormalQMC	基于 QMC 的多元正态随机数

Helpers（辅助函数）

函数	说明
discrepancy	计算点集在超立方体中的差异度（均匀性指标）
geometric_discrepancy	基于几何性质（最小距离 / 最小生成树）的均匀性指标
update_discrepancy	增量式更新 centered discrepancy
scale	在单位超立方体与任意边界之间做缩放变换

三、点集质量度量：discrepancy 三件套

3.1 discrepancy：四种经典差异度

函数签名（见 scipy/stats/_qmc.py 第 197-202 行）：

discrepancy(sample, *, iterative=False, method="CD", workers=1)

概念：discrepancy 是衡量超立方体内点集填充均匀性的准则，量化的是“超立方体上的连续均匀分布”与“n 个不同采样点上的离散均匀分布”之间的距离。值越低，参数空间覆盖越好。对超立方体的一组子集，discrepancy 等于落在某个子集中的样本点比例与该子集体积之差（部分定义取各子集上的均方根差而非最大值）。

模块源码指出，合理的均匀性度量应满足：因子/运行置换下不变、坐标旋转下不变、能度量全空间与低维投影均匀性、有合理几何意义、易计算、满足 Koksma-Hlawka 型不等式、与试验设计其他准则一致。

四种可用方法：

• CD（Centered Discrepancy，居中差异）：子空间涉及超立方体的一个角；默认值。实现对应 Zhou 等论文的式 9 右侧（不取平方根）。
• WD（Wrap-around Discrepancy，环绕差异）：子空间可环绕边界；对应式 10。
• MD（Mixture Discrepancy，混合差异）：CD/WD 的混合，覆盖更多准则；对应式 18。
• L2-star：类似 CD 但对旋转不变；按 Warnock 的式 10 实现，取平方根。

在源码实现中，四种方法分别映射到 _qmc_cy 中的 Cython 封装（第 321-326 行）：

methods = {
    "CD": _cy_wrapper_centered_discrepancy,
    "WD": _cy_wrapper_wrap_around_discrepancy,
    "MD": _cy_wrapper_mixture_discrepancy,
    "L2-star": _cy_wrapper_l2_star_discrepancy,
}

iterative=True 的用途：按“已有 n 个样本、假装再多 1 个”的方式计算差异度。当你想给点集新增一个点、且需要比较多个候选点时，可以先对初始样本算一次 iterative 差异度，再对每个候选点调用 update_discrepancy 增量更新——这比逐候选重算整组差异度快得多。

workers：并行处理的 worker 数，传 -1 表示使用所有 CPU 线程，默认 1。相关测试见 scipy/stats/tests/test_qmc.py 的 test_discrepancy_parallel。

官方示例（计算 6 个点在边界 [0.5, 6.5] 内数据的差异度）：

import numpy as np
from scipy.stats import qmc

space = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
l_bounds = [0.5, 0.5]
u_bounds = [6.5, 6.5]
space = qmc.scale(space, l_bounds, u_bounds, reverse=True)
print(space)
# array([[0.08333333, 0.41666667],
#        [0.25      , 0.91666667],
#        [0.41666667, 0.25      ],
#        [0.58333333, 0.75      ],
#        [0.75      , 0.08333333],
#        [0.91666667, 0.58333333]])
print(qmc.discrepancy(space))
# 0.008142039609053464

3.2 geometric_discrepancy：几何视角的均匀性

函数签名：

geometric_discrepancy(sample, method="mindist", metric="euclidean")

它基于点集的几何性质度量样本质量：method="mindist"（默认）取任意两点间的最小距离，method="mst" 取最小生成树（minimum spanning tree）的平均边长。metric 为距离度量，可用 scipy.spatial.distance.pdist 支持的任何度量。

注意它与 discrepancy 的方向相反：这里值越大表示覆盖越好（点间距离越大越均匀），而 discrepancy 是值越小越好。文档同时提醒：比较不同采样策略时，样本量必须保持不变；mst 方法目前只计算平均边长（最小生成树还可导出边长标准差，均值更高、标准差更低更优，见 Franco 等的论文讨论）。

从实现看，mindist 返回的是 distances[distances.nonzero()] 的最小值，遇到重复点会发出 Sample contains duplicate points. 警告；mst 借助 scipy.sparse.csgraph.minimum_spanning_tree 与 scipy.spatial.distance.squareform 完成（scipy/stats/_qmc.py 第 447-455 行）。

rng = np.random.default_rng(191468432622931918890291693003068437394)
sample = qmc.LatinHypercube(d=2, rng=rng).random(50)
print(qmc.geometric_discrepancy(sample))            # 0.03708161435687876 (mindist)
print(qmc.geometric_discrepancy(sample, method='mst'))  # 0.1105149978798376

3.3 update_discrepancy：增量更新

update_discrepancy(x_new, sample, initial_disc)

给定初始样本 sample (n, d) 及其 centered discrepancy initial_disc，加入新点 x_new (1, d) 后返回更新后的差异度。典型用法（与 discrepancy(..., iterative=True) 配合）：

disc_init = qmc.discrepancy(space[:-1], iterative=True)
print(disc_init)                              # 0.04769081147119336
print(qmc.update_discrepancy(space[-1], space[:-1], disc_init))
# 0.008142039609053513  （与直接计算整组一致）

实现直接调用 Cython 封装 _cy_wrapper_update_discrepancy（scipy/stats/_qmc.py 第 519 行）。值得留意的是，同一文件中的 _perturb_discrepancy 实现了 Jin 等论文的“基本扰动”方案（交换两行某列坐标 sample[i1, k] <-> sample[i2, k] 后 O(1) 更新差异度），拉丁超立方优化正是依赖这类技巧实现的。

四、核心引擎详解

4.1 QMCEngine：统一接口

QMCEngine 是所有引擎的抽象基类，定义了统一的采样与状态管理接口（scipy/stats/_qmc.py 第 799 行起）：

• random(n=1, *, workers=1)：生成 n 行样本；
• reset()：重置引擎到初始状态，返回 self，便于复现；
• fast_forward(n)：跳过 n 个点后返回 self；
• integers(l_bounds, *, u_bounds=None, n=1, endpoint=False, workers=1)：生成整数样本，用于离散问题。

你还可以继承 QMCEngine 实现自定义引擎（测试见 test_qmc.py 的 test_subclassing_QMCEngine），只需实现 _random 并调用基类的输入校验即可。

4.2 Sobol：最广泛使用的低差异序列

构造函数：

Sobol(d, *, scramble=True, bits=None, rng=None, optimization=None)

• d：维度，方向数来自 Joe & Kuo 的构造，最多支持 21201 维；
• scramble=True：是否打乱（scrambled Sobol' 属于 RQMC，具备前文所述的方差与稳健性优势）；
• bits：生成点的位宽，可为 32/64（默认自动选择），影响引擎的状态空间；
• optimization：可选的优化后处理（见下文“优化”小节）；
• rng：随机数生成器（打乱与优化需要）。

特殊方法 random_base2(m)：一次生成恰好 2^m 个点。由于 Sobol' 的特殊样本量是 2 的幂，random_base2 能保证样本的完整结构性质（test 中 test_random_base2 验证了这一点）；需要约 2^m 个点时用它优于 random(n)。其余用法与 QMCEngine 一致（reset、fast_forward 均已覆写以同步打乱状态）。

from scipy.stats import qmc

sampler = qmc.Sobol(d=2, scramble=False)
sample = sampler.random_base2(m=3)   # 恰好 8 个点
print(sample.shape)                  # (8, 2)
sampler.reset()                      # 重新开始

4.3 Halton：按素数基构造

Halton(d, *, scramble=True, optimization=None, rng=None)

Halton 序列用前 d 个素数作为基分别构造 van der Corput 序列再拼接而成。底层单维序列由 van_der_corput 函数生成（scipy/stats/_qmc.py 第 730 行），其签名：

van_der_corput(n, base=2, *, start_index=0, scramble=False,
               permutations=None, rng=None, workers=1)

permutations 可传入预生成的基置换（用 _van_der_corput_permutations(base, rng=...) 生成）以实现可复现的打乱；start_index 允许从序列中某个位置开始。

如前文所述，Halton 序列前几个维度的均分性质远好于后面的维度，且没有特殊样本量约束，通常精度不及 Sobol'。

4.4 LatinHypercube：试验设计的主流选择

LatinHypercube(d, *, scramble=True, strength=1, optimization=None, rng=None)

拉丁超立方采样（LHS）将每个维度划分为 n 个等宽区间，保证每行每列恰好一个点。关键参数：

• scramble：若为 False，点集完全确定（每维取区间中点）；若为 True（默认），在保持 LHS 结构的前提下随机打乱；
• strength：正交性强度，strength=1 为经典 LHS；strength=2 为 OA-LHS（基于正交数组的 LHS，实现为 _random_oa_lhs），在低维投影上更均匀；更高强度仅当 n 为 4 的幂时可用（源码中 _random_oa_lhs(n=4) 暗示了这一约束，相关测试为 test_sample_stratified）。

LHS 与 Sobol' 的一个核心区别：LHS 对任意 n 都适用（不存在“特殊样本量”问题），因此当点数不是 2 的幂时它通常是更自然的选择。

4.5 PoissonDisk：保证最小间距的点集

PoissonDisk(d, *, radius=0.05, hypersphere="volume", ncandidates=30,
            optimization=None, rng=None, l_bounds=None, u_bounds=None)

泊松盘采样保证任意两点间距离不小于 radius，适合需要“空间填充 + 间距约束”的场景：

• radius：点间最小距离（默认 0.05）；
• hypersphere："volume"（默认）或 "surface"，决定候选点在以已有点为中心的超球内是“体积采样”还是“表面采样”；
• ncandidates：每步尝试的候选点数量（默认 30）；
• l_bounds / u_bounds：与 scale 类似，可指定采样边界（默认单位超立方体）。

实现内部维护网格池 _initialize_grid_pool 加速邻域查询，并提供 fill_space() 方法持续填充直到空间填满（测试见 test_fill_space）。

4.6 MultivariateNormalQMC 与 MultinomialQMC：向真实分布前进

MultivariateNormalQMC：用 QMC 点生成多元正态随机数。

MultivariateNormalQMC(mean, cov=None, *, cov_root=None, inv_transform=True,
                      engine=None, rng=None)

• mean：均值；cov：协方差（或用 cov_root 直接给 Cholesky 因子）；
• inv_transform=True：用逆变换法（对标准正态做均匀→正态分位数变换），这是默认且推荐的模式；设为 False 则用 Box-Muller 类方法；
• engine：底层 QMC 引擎，默认是打乱的 Sobol；可传入任意 QMCEngine（测试 test_other_engine 验证了这一点）。

engine = qmc.MultivariateNormalQMC(mean=[0.1, 0.5], cov=[[1, 0.3], [0.3, 1]])
samples = engine.random(n=100)      # (100, 2) 的多元正态样本

MultinomialQMC：用 QMC 点生成多项分布随机数。

MultinomialQMC(pvals, n_trials, *, engine=None, rng=None)

pvals 为各类概率（和需为 1），n_trials 为试验次数；内部将 QMC 均匀点按累积概率分箱映射为类别计数（Cython 辅助 _fill_p_cumulative、_categorize），保证行和为 n_trials。测试 test_MultinomialBasicDraw、test_MultinomialDistribution 覆盖其正确性。

五、scale：单位超立方体与真实边界的桥梁

QMC 引擎输出恒位于 [0, 1]^d，真实问题几乎总是需要其他边界。scale 函数完成这一变换（scipy/stats/_qmc.py 第 84-164 行）：

scale(sample, l_bounds, u_bounds, *, reverse=False)

正变换（默认）把 [0, 1) 映射到 [a, b)，公式为 (b - a) * sample + a；reverse=True 则把位于边界内的点映射回单位超立方体。实现会校验：正变换要求样本在单位超立方体内，反变换要求样本在边界内，且样本必须是 2D 数组。

l_bounds = [-2, 0]
u_bounds = [6, 5]
sample = [[0.5, 0.75],
          [0.5, 0.5],
          [0.75, 0.25]]
sample_scaled = qmc.scale(sample, l_bounds, u_bounds)
print(sample_scaled)
# array([[2.  , 3.75],
#        [2.  , 2.5 ],
#        [4.  , 1.25]])
back = qmc.scale(sample_scaled, l_bounds, u_bounds, reverse=True)
print(back)
# array([[0.5 , 0.75],
#        [0.5 , 0.5 ],
#        [0.75, 0.25]])

六、优化后处理：random-cd 与 lloyd

LatinHypercube、Sobol、Halton、PoissonDisk 均支持 optimization 参数（None、"random-cd" 或 "lloyd"），由模块内部 _select_optimizer 分发（scipy/stats/_qmc.py 第 2585 行）：

• "random-cd"：随机化 centered discrepancy 优化。随机选取两行做“基本扰动”（交换第 k 列坐标），若扰动降低差异度则接受，重复直到满足停止条件。它借助 _perturb_discrepancy 的 O(1) 增量更新，避免重算整个点集的差异度——这正是 Jin 等论文算法的 SciPy 实现。LatinHypercube 文档说明：使用 random-cd 时每次调用 random 生成的点数需大于 1 才有意义。
• "lloyd"：Lloyd 质心 Voronoi 迭代（CVT）。每步把每个点移动到其 Voronoi 胞腔的质心（_lloyd_iteration 借助 scipy.spatial.Voronoi），_lloyd_centroidal_voronoi_tessellation 默认 tol=1e-5、maxiter=10。测试 test_optimizers、test_lloyd、test_lloyd_non_mutating 覆盖其行为。

两者都只在非打乱或打乱后对点集做后处理；optimization 与 scramble=True 并用时，优化作用于每次生成的整批点（测试 test_optimizers 验证了不同优化对差异度/几何差异度指标的影响）。

七、快速上手：一个完整的 QMC 数值积分示例

综合以上内容，一个典型的 QMC 工作流是：选引擎 → 生成单位立方体点 → 缩放/变换 → 计算被积函数均值 → 用 discrepancy 或 RQMC 重复评估精度。这里演示一个经典的例子——估计三维单位球内指示函数的积分（真实值 4π/3 / 8 = π/6 ≈ 0.5236，该指示函数不满足 Hardy-Krause 光滑性要求，但足以对比 MC 与 QMC 的方差）：

import numpy as np
from scipy.stats import qmc

def mc_estimate(n, d=3, seed=42):
    rng = np.random.default_rng(seed)
    x = rng.random((n, d))
    return np.mean(np.sum(x**2, axis=1) <= 1.0)

def qmc_estimate(n, d=3, scramble=True):
    engine = qmc.Sobol(d=d, scramble=scramble)
    x = engine.random_base2(int(np.log2(n))) if n & (n-1) == 0 else engine.random(n)
    return np.mean(np.sum(x**2, axis=1) <= 1.0)

for n in [64, 256, 1024]:
    print(n,
          "MC :", mc_estimate(n),
          "QMC:", qmc_estimate(n))

注意：为保持 Sobol 的特殊样本量性质，上例在 n 为 2 的幂时使用 random_base2。若 n 不是 2 的幂，则应改用 LatinHypercube（任意 n 都适用）或接受 Sobol 点集的轻微结构退化（这正是模块文档强调“改变样本量会降低性能”的体现）。

八、进一步阅读与仓库定位

• 模块入口与公开 API：from scipy.stats import qmc，__all__ 定义于 scipy/stats/_qmc.py 第 44-46 行；
• 模块级数学背景文档：scipy/stats/qmc.py 顶部 docstring（含 27 篇参考文献，覆盖 Owen、Niederreiter、Dick/Kuo/Sloan、Joe & Kuo 方向数等）；
• 参考文档页面：doc/source/reference/stats.qmc.rst；
• 单元测试：scipy/stats/tests/test_qmc.py 覆盖 scale、四种 discrepancy、update_discrepancy、van_der_corput、各引擎的采样/继续/重置/快进/边界、Sobol 64 位、MultivariateNormalQMC 的 Shapiro 正态性检验、Lloyd 优化等；
• Cython 加速内核位于 scipy/stats/_qmc_cy.pyx（差异度、van der Corput、多项分布采样等热路径）。

适用前提与限制：scipy.stats.qmc 的引擎均基于 NumPy 数组在 CPU 上运行；Sobol 维度上限 21201、random_base2 要求 2 的幂、strength=2 的 OA-LHS 要求特定 n 等约束来自当前仓库源码，使用时需按版本与官方文档核对。本文所有示例均可直接在装有 SciPy 的环境中运行验证，无需修改仓库任何文件。