拟蒙特卡罗#

在讨论拟蒙特卡罗 (QMC) 之前，先简要介绍一下蒙特卡罗 (MC)。MC 方法（或称 MC 实验）是一类广泛的计算算法，依靠重复随机采样来获得数值结果。其核心理念是利用随机性来解决原则上可能是确定性的问题。它们常用于物理和数学问题，在其他方法难以实现或无法使用时最为有用。MC 方法主要用于三类问题：优化、数值积分以及从概率分布中生成抽取样本。

生成具有特定属性的随机数是一个比听起来更复杂的问题。简单的 MC 方法旨在采样独立同分布 (IID) 的点。但是，生成多组随机点可能会产生截然不同的结果。

在上图的两种情况下，点的生成都是随机的，不了解之前抽取的任何点。显而易见，空间的某些区域未被探索——这可能会在模拟中导致问题，因为特定的一组点可能会触发完全不同的行为。

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})\)。

虽然收敛性得到了保证，但从业者往往希望拥有一个更具确定性的探索过程。在普通 MC 中，可以使用种子来实现可重复的过程。但固定种子会破坏收敛特性：给定的种子可能适用于某一类问题，但在另一类问题上失效。

为了以确定性的方式遍历空间，通常的做法是使用跨越所有参数维度的正则网格，也称为饱和设计 (saturated design)。假设一个单位超立方体，所有边界范围都在 0 到 1 之间。现在，如果点之间的距离为 0.1，填满单位区间所需的点数将是 10 个。在 2 维超立方体中，同样的间距需要 100 个点，而在 3 维中则需要 1,000 个点。随着维度增加，填满空间所需的实验次数随空间维度的增加而呈指数级增长。这种指数级增长被称为“维度灾难”。

>>> 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 方法应运而生。它们是确定性的，具有良好的空间覆盖性，且其中一些方法可以持续采样并保持良好的特性。与 MC 方法的主要区别在于，这些点不是独立同分布的，而是“知道”之前的点。因此，某些方法也被称为序列 (sequences)。

该图展示了 2 组 256 个点。左侧的设计是普通的 MC，而右侧的设计是使用 Sobol’ 方法的 QMC 设计。我们可以清楚地看到 QMC 版本更加均匀。这些点在边界附近的采样效果更好，且聚集簇或间隙更少。

评估均匀性的一种方法是使用一种称为“偏差 (discrepancy)”的度量。在这里，Sobol’ 点的偏差优于粗略的 MC。

回到均值的计算，QMC 方法在误差收敛速率上也表现更好。对于此函数，它们可以达到 \(O(n^{-1})\)，在非常平滑的函数上甚至会有更好的速率。此图显示 Sobol’ 方法的速率为 \(O(n^{-1})\)。

更多数学细节请参考 scipy.stats.qmc 的文档。

计算偏差#

让我们考虑两组点。从下图可以清楚地看到，左侧的设计比右侧的设计覆盖了更多的空间。这可以使用 scipy.stats.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 引擎#

已经实现了几种 QMC 采样器/引擎。这里我们看看两种最常用的 QMC 方法：scipy.stats.qmc.Sobol 和 scipy.stats.qmc.Halton 序列。

"""Sobol' and Halton sequences."""
from scipy.stats import qmc
import numpy as np

import matplotlib.pyplot as plt


rng = np.random.default_rng()

n_sample = 256
dim = 2

sample = {}

# Sobol'
engine = qmc.Sobol(d=dim, seed=rng)
sample["Sobol'"] = engine.random(n_sample)

# Halton
engine = qmc.Halton(d=dim, seed=rng)
sample["Halton"] = engine.random(n_sample)

fig, axs = plt.subplots(1, 2, figsize=(8, 4))

for i, kind in enumerate(sample):
    axs[i].scatter(sample[kind][:, 0], sample[kind][:, 1])

    axs[i].set_aspect('equal')
    axs[i].set_xlabel(r'$x_1$')
    axs[i].set_ylabel(r'$x_2$')
    axs[i].set_title(f'{kind}—$C^2 = ${qmc.discrepancy(sample[kind]):.2}')

plt.tight_layout()
plt.show()

警告

QMC 方法需要特别注意，用户必须阅读文档以避免常见陷阱。例如，Sobol’ 要求点数必须是 2 的幂。此外，抽稀 (thinning)、烧入 (burning) 或其他点选择操作可能会破坏序列的特性，导致结果点集并不优于 MC。

QMC 引擎具有状态感知能力。这意味着你可以继续序列、跳过某些点或重置它。让我们从 scipy.stats.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]])

注意

默认情况下，scipy.stats.qmc.Sobol 和 scipy.stats.qmc.Halton 都是经过加乱 (scrambled) 处理的。这样收敛特性更好，并能防止在高维空间中出现条纹或明显的点阵模式。在实际应用中，没有理由不使用加乱版本。

构建 QMC 引擎，即继承 `QMCEngine`#

要创建你自己的 scipy.stats.qmc.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]])

QMC 使用指南#

QMC 是有规则的！请务必阅读文档，否则你可能无法获得优于 MC 的益处。
如果你需要恰好 \(2^m\) 个点，请使用 scipy.stats.qmc.Sobol。
scipy.stats.qmc.Halton 允许采样或跳过任意数量的点。这是以收敛速率慢于 Sobol’ 为代价的。
切勿删除序列的前几个点。这会破坏其特性。
加乱 (Scrambling) 总是更好的选择。
如果你使用基于拉丁超立方采样 (LHS) 的方法，则无法在不丢失 LHS 特性的情况下添加点。（虽然有一些方法可以做到这一点，但目前尚未实现。）

拟蒙特卡罗#

计算偏差#

使用 QMC 引擎#

构建 QMC 引擎，即继承 QMCEngine#

QMC 使用指南#

构建 QMC 引擎，即继承 `QMCEngine`#