Python SciPy qmc.Sobol用法及代码示例

本文简要介绍 python 语言中 scipy.stats.qmc.Sobol 的用法。

用法:

class  scipy.stats.qmc.Sobol(d, *, scramble=True, bits=None, seed=None, optimization=None)#

用于生成(加扰)Sobol 序列的引擎。

Sobol 的序列是 low-discrepancy、quasi-random 数字。可以使用两种方法绘制点：

• random_base2 ：安全地绘制\(n=2^m\) 点。这种方法保证了序列的平衡特性。

• random ：从序列中绘制任意数量的点。请参阅下面的警告。

参数 ：：

d： int

序列的维度。最大维数为 21201。

scramble： 布尔型，可选

如果为 True，则使用 LMS+shift 加扰。否则，不进行加扰。默认为真。

bits： 整数，可选

生成器的位数。控制可以生成的最大点数，即 2**bits 。最大值为 64。它与返回类型不对应，返回类型始终为np.float64，以防止点重复。默认值为 None，为了向后兼容，对应于 30。

optimization： {无，“random-cd”，“lloyd”}，可选

是否使用优化方案来提高采样后的质量。请注意，这是一个后处理步骤，不能保证样品的所有属性都得到保留。默认为“无”。

• random-cd：坐标的随机排列以降低中心差异。基于中心差异的最佳样本不断更新。与使用其他差异度量相比，基于中心差异的采样对 2D 和 3D 子投影显示出更好的 space-filling 鲁棒性。

• lloyd：使用修改后的Lloyd-Max 算法扰动样本。该过程收敛到等间隔的样本。

seed： {无，整数， numpy.random.Generator }，可选

如果种子是 int 或 None，一个新的numpy.random.Generator是使用创建的np.random.default_rng(seed).如果种子已经是一个Generator实例，然后使用提供的实例。

注意：

Sobol 序列 [1] 在 \([0,1)^{d}\) 中提供了 \(n=2^m\) 低差异点。对它们进行置乱[3]使它们适用于奇异被积函数，提供了一种误差估计方法，并且可以提高它们的收敛速度。所实施的加扰策略是(左)线性矩阵加扰(LMS)，然后是数字随机移位(LMS+移位)[2]。

索博尔序列有许多版本，具体取决于它们的“方向号”。该代码使用[4]中的方向号。因此，最大维度数为 21201。方向数已使用搜索标准 6 预先计算，并且可以在 https://web.maths.unsw.edu.au/~fkuo/sobol/ 处检索。

警告

Sobol 序列是一个求积规则，如果使用的样本量不是 2 的幂，或者跳过第一个点，或者对序列进行细化 [5]，它们就会失去平衡特性。

如果 \(n=2^m\) 点数不够，则应将 \(2^M\) 点数用于 \(M>m\) 。加扰时，独立重复的数量 R 不必是 2 的幂。

Sobol' 序列被生成到一定数量\(B\)位。后\(2^B\)点已经生成，序列将重复。因此，出现了错误。位数可以通过参数控制位.

参考：

[1]

I. M. Sobol，“立方体中点的分布和积分的准确评估。” ZH。维奇斯尔。填充。我在。物理学，7：784-802，1967。

[2]

J. Matousek，“关于锚定框的 L2 差异”。 J. of Complexity 14, 527-556, 1998。

[3]

Art B. Owen，“加扰 Sobol 和Niederreiter-Xing 点”。复杂性杂志，14(4):466-489，1998 年 12 月。

[4]

S. Joe 和 F. Y. Kuo，“用更好的二维投影构建 sobol 序列”。 SIAM 科学计算杂志，30(5):2635-2654，2008。

[5]

Art B. Owen，“关于放弃第一个 Sobol 点。” arXiv:2008.08051，2020 年。

例子：

从 Sobol 的低差异序列中生成样本。

>>> from scipy.stats import qmc
>>> sampler = qmc.Sobol(d=2, scramble=False)
>>> sample = sampler.random_base2(m=3)
>>> sample
array([[0.   , 0.   ],
       [0.5  , 0.5  ],
       [0.75 , 0.25 ],
       [0.25 , 0.75 ],
       [0.375, 0.375],
       [0.875, 0.875],
       [0.625, 0.125],
       [0.125, 0.625]])

使用差异标准计算样本的质量。

>>> qmc.discrepancy(sample)
0.013882107204860938

要继续现有设计，可以通过再次调用 random_base2 来获得加分。或者，您可以跳过一些要点，例如：

>>> _ = sampler.reset()
>>> _ = sampler.fast_forward(4)
>>> sample_continued = sampler.random_base2(m=2)
>>> sample_continued
array([[0.375, 0.375],
       [0.875, 0.875],
       [0.625, 0.125],
       [0.125, 0.625]])

最后，样本可以缩放到边界。

>>> l_bounds = [0, 2]
>>> u_bounds = [10, 5]
>>> qmc.scale(sample_continued, l_bounds, u_bounds)
array([[3.75 , 3.125],
       [8.75 , 4.625],
       [6.25 , 2.375],
       [1.25 , 3.875]])