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

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

用法:

class  scipy.stats.qmc.LatinHypercube(d, *, scramble=True, strength=1, optimization=None, seed=None)#

拉丁超立方抽样 (LHS)。

拉丁超立方体样本 [1] 在 \([0,1)^{d}\) 中生成 \(n\) 点。每个单变量边际分布都是分层的，在 \([j/n, (j+1)/n)\) 中为 \(j=0,1,...,n-1\) 放置一个点。它们在 \(n << d\) 时仍然适用。

参数 ：：

d： int

参数空间的维度。

scramble： 布尔型，可选

当为 False 时，将样本集中在多维网格的单元格内。否则，样本会随机放置在网格的单元格内。

注意

环境scramble=False不确保确定性输出。为此，请使用种子范围。

默认为 True。

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

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

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

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

strength： {1, 2}，可选

LHS 的强度。 strength=1 产生一个普通的 LHS，而strength=2 产生一个基于正交阵列的强度为 2 [7]、[8] 的 LHS。在这种情况下，只能对n=p**2 点进行采样，而p 是一个素数。它还约束 d <= p + 1 。默认值为 1。

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

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

注意：

当 LHS 用于对函数 \(f\) 对 \(n\) 进行积分时，LHS 对于几乎可加的被积数非常有效 [2]。对于 \(n\) 点的 LHS，积分方差始终低于 \(n-1\) 点上的普通 MC [3]。 LHS 在积分的均值和方差上存在中心极限定理 [4]，但由于随机化，优化的 LHS 不一定有中心极限定理。

如果在 \(A\) 的每个 n-row-by-t 列子矩阵中，\(A\) 被称为强度 \(t\) 的正交数组：所有 \(p^t\) 可能的不同行出现相同的次数。 \(A\) 的元素位于集合 \(\{0, 1, ..., p-1\}\) 中，也称为符号。 \(p\) 必须是素数的约束是为了允许模运算。增加强度会增加样本 sub-projections 的一些对称性。强度为 2 时，样本沿 2D sub-projections 的对角线对称。这可能是不希望的，但另一方面，样品分散得到改善。

强度 1(普通 LHS)比强度 0 (MC) 具有优势，强度 2 是比强度 1 有用的增量。强度 3 是一个较小的增量，并且像 Sobol'、Halton 这样的加扰 QMC 性能更高 [7]。

为了创建强度为 2 的 LHS，正交数组 \(A\) 通过将符号集的随机双射映射应用于自身来随机化。例如，在第 0 列中，所有 0 都可能变为 2；在第 1 列中，所有 0 都可能变为 1，依此类推。然后，对于每一列 \(i\) 和符号 \(j\) ，我们将大小为 \(p\) 的普通一维 LHS 添加到 \(A^i = j\) 的子数组中。结果矩阵最终除以 \(p\) 。

参考：

[1]

Mckay 等人，“在计算机代码输出分析中选择输入变量值的三种方法的比较”。技术计量学，1979 年。

[2]

M. Stein，“使用拉丁超立方抽样的模拟的大样本属性”。技术计量学 29，没有。 2：143-151，1987。

[3]

A. B. Owen，“加扰网络正交的蒙特卡洛方差”。 SIAM 数值分析杂志 34，没有。 5: 1884-1910, 1997

[4]

罗，W.-L。 “关于拉丁超立方抽样。”统计年鉴24，没有。 5：2058-2080，1996。

[5]

方等人。 “计算机实验的设计和建模”。计算机科学与数据分析系列，2006 年。

[6]

Damblin 等人，“空间填充设计的数值研究：拉丁超立方体样本和子投影属性的优化”。模拟杂志，2013 年。

[7] (1,2)

A. B. Owen，“用于计算机实验、集成和可视化的正交阵列”。中央统计局，1992 年。

[8]

B. Tang，“基于正交阵列的拉丁超立方体”。美国统计协会杂志，1993 年。

[9]

苏珊·K·西霍姆等人。 “拉丁超立方抽样和蒙特卡罗流行病模型的敏感性分析”。 Int 生物医学计算杂志，23(1-2), 97-112，DOI:10.1016/0020-7101(88)90067-0，1988。

例子：

在[9]中，使用拉丁超立方采样策略对参数空间进行采样，以研究流行病模型每个参数的重要性。这种分析也称为敏感性分析。

由于问题的维数很高 (6)，因此覆盖该空间的计算成本很高。当数值实验成本高昂时，QMC 可以实现使用网格可能无法进行的分析。

模型的六个参数代表患病概率、退出概率和四个接触概率。作者假设所有参数均匀分布并生成 50 个样本。

使用scipy.stats.qmc.LatinHypercube复制协议，第一步是在单位超立方体中创建样本：

>>> from scipy.stats import qmc
>>> sampler = qmc.LatinHypercube(d=6)
>>> sample = sampler.random(n=50)

然后可以将样本缩放到适当的范围：

>>> l_bounds = [0.000125, 0.01, 0.0025, 0.05, 0.47, 0.7]
>>> u_bounds = [0.000375, 0.03, 0.0075, 0.15, 0.87, 0.9]
>>> sample_scaled = qmc.scale(sample, l_bounds, u_bounds)

使用这样的样本运行模型 50 次，并构建了多项式响应面。这使得作者能够研究每个参数在每个其他参数的可能性范围内的相对重要性。在这项计算机实验中，与网格采样相比，他们将响应曲面上的误差保持在 2% 以下所需的样本数量减少了 14 倍。

以下是其他示例，展示了构建具有更好空间覆盖范围的 LHS 的替代方法。

使用基本 LHS 作为基线。

>>> sampler = qmc.LatinHypercube(d=2)
>>> sample = sampler.random(n=5)
>>> qmc.discrepancy(sample)
0.0196...  # random

使用优化关键字参数以较高的计算成本产生具有较小差异的 LHS。

>>> sampler = qmc.LatinHypercube(d=2, optimization="random-cd")
>>> sample = sampler.random(n=5)
>>> qmc.discrepancy(sample)
0.0176...  # random

使用强度关键字参数来生成强度为 2 的基于正交数组的 LHS。在这种情况下，样本点的数量必须是素数的平方。

>>> sampler = qmc.LatinHypercube(d=2, strength=2)
>>> sample = sampler.random(n=9)
>>> qmc.discrepancy(sample)
0.00526...  # random

可以组合选项以产生优化的基于中心正交阵列的 LHS。优化后的结果不能保证强度为2。