scipy.stats.qmc.

拉丁超立方#

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整数

参数空间的维度。

打乱bool, 可选

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

备注

设置 scramble=False 并不能确保输出是确定性的。为此,请使用 seed 参数。

默认值为 True。

Added in version 1.10.0.

优化{None, “random-cd”, “lloyd”}, 可选

是否使用优化方案来提高采样后的质量。请注意,这是一个后处理步骤,不保证样本的所有属性都将被保留。默认值为 None。

  • random-cd: 坐标随机排列以降低中心差异。基于中心差异的最佳样本不断更新。与使用其他差异度量相比,基于中心差异的采样在2D和3D子投影方面表现出更好的空间填充鲁棒性。

  • lloyd: 使用改进的Lloyd-Max算法扰动样本。该过程收敛于等间距的样本。

Added in version 1.8.0.

在 1.10.0 版本发生变更: 添加 lloyd

强度{1, 2}, 可选

LHS 的强度。strength=1 生成一个普通的 LHS,而 strength=2 生成一个基于强度为 2 的正交数组的 LHS [7][8]。在这种情况下,只能采样 n=p**2 个点,其中 p 是一个质数。它还限制了 d <= p + 1。默认值为 1。

Added in version 1.8.0.

种子 : {None, int, numpy.random.Generator}, 可选{None, int,}

如果 seed 是 int 或 None,将使用 np.random.default_rng(seed) 创建一个新的 numpy.random.Generator。如果 seed 已经是 Generator 实例,则使用提供的实例。

方法

fast_forward(n)

将序列快进 n 个位置。

integers(l_bounds, *[, u_bounds, n, ...])

l_bounds`(包含)到 `u_bounds`(不包含)之间抽取 `n 个整数,或者如果 endpoint=True,则从 `l_bounds`(包含)到 `u_bounds`(包含)之间抽取。

random([n, workers])

在半开区间 [0, 1) 中绘制 n

reset()

将引擎重置为基础状态。

注释

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

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

强度 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, “Large sample properties of simulations using Latin hypercube sampling.” Technometrics 29, no. 2: 143-151, 1987.

[3]

A. B. Owen, “Monte Carlo variance of scrambled net quadrature.” SIAM Journal on Numerical Analysis 34, no. 5: 1884-1910, 1997

[4]

Loh, W.-L. “关于拉丁超立方抽样。” 《统计年鉴》24, 第5期: 2058-2080, 1996.

[5]

Fang 等人。《计算机实验的设计与建模》。计算机科学与数据分析系列,2006年。

[6]

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

[7] (1,2)

A. B. Owen , “Orthogonal arrays for computer experiments, integration and visualization.” Statistica Sinica, 1992.

[8]

B. Tang, “Orthogonal Array-Based Latin Hypercubes.” Journal of the American Statistical Association, 1993.

[9]

Susan K. Seaholm 等人。“拉丁超立方抽样与蒙特卡罗流行病模型的敏感性分析”。《国际生物医学计算杂志》,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

使用 optimization 关键字参数以更高的计算成本生成具有较低差异性的 LHS。

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

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

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

选项可以组合以生成基于优化的中心正交数组的LHS。优化后,结果不能保证具有强度2。