scipy.stats.qmc.

LatinHypercube#

class scipy.stats.qmc.LatinHypercube(d, *, scramble=True, strength=1, optimization=None, rng=None)[source]#

拉丁超立方抽样 (LHS)。

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

参数:
dint

参数空间的维度。

scramblebool,可选

如果为 False,则将样本居中放置在多维网格的单元格中。否则,样本会随机放置在网格的单元格中。

注意

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

默认为 True。

1.10.0 版本中新增。

optimization{None, “random-cd”, “lloyd”}, optional

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

  • random-cd:随机置换坐标以降低中心偏差。基于中心偏差的最佳样本会不断更新。与其他偏差测量方法相比,基于中心偏差的采样在 2D 和 3D 子投影中显示出更好的空间填充鲁棒性。

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

1.8.0 版本中新增。

在 1.10.0 版本中更改: 添加 lloyd

strength{1, 2}, optional

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

1.8.0 版本中新增。

rngnumpy.random.Generator, optional

伪随机数生成器状态。当 rng 为 None 时,将使用来自操作系统的熵创建一个新的 numpy.random.Generatornumpy.random.Generator 以外的类型将传递给 numpy.random.default_rng 以实例化一个 Generator

在 1.15.0 版本中更改: 作为从使用 numpy.random.RandomState 过渡到 numpy.random.GeneratorSPEC-007 的一部分,此关键字已从 seed 更改为 rng。在过渡期内,两个关键字将继续有效,但一次只能指定一个关键字。过渡期结束后,使用 seed 关键字的函数调用将发出警告。经过弃用期后,将删除 seed 关键字。

方法

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 用于在 \(n\) 上积分函数 \(f\) 时,LHS 在几乎是可加性的被积函数上非常有效 [2]。对于 \(n\) 个点的 LHS,积分的方差始终低于 \(n-1\) 个点上的普通 MC [3]。LHS 在积分的均值和方差上有一个中心极限定理 [4],但由于随机化,对于优化的 LHS 不一定成立。

如果 \(A\) 的每个 n 行 t 列子矩阵中:所有 \(p^t\) 个可能的不同行出现次数相同,则 \(A\) 被称为强度为 \(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 等人,“比较三种在计算机代码输出分析中选择输入变量值的方法”。Technometrics,1979。

[2]

M. Stein,“使用拉丁超立方采样的模拟的大样本属性”。Technometrics 29,no. 2: 143-151, 1987。

[3]

A. B. Owen,“加扰网格正交的蒙特卡罗方差”。SIAM Journal on Numerical Analysis 34, no. 5: 1884-1910, 1997

[4]

Loh, W.-L.“关于拉丁超立方采样”。统计年鉴 24, no. 5: 2058-2080, 1996。

[5]

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

[6]

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

[7] (1,2)

A. B. Owen,“用于计算机实验、积分和可视化的正交阵列”。Statistica Sinica,1992 年。

[8]

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

[9]

Seaholm, Susan K. 等人 (1988)。拉丁超立方抽样和蒙特卡罗流行病模型的灵敏度分析。Int J Biomed Comput, 23(1-2), 97-112。 DOI:10.1016/0020-7101(88)90067-0

示例

从拉丁超立方生成器生成样本。

>>> from scipy.stats import qmc
>>> sampler = qmc.LatinHypercube(d=2)
>>> sample = sampler.random(n=5)
>>> sample
array([[0.1545328 , 0.53664833], # random
        [0.84052691, 0.06474907],
        [0.52177809, 0.93343721],
        [0.68033825, 0.36265316],
        [0.26544879, 0.61163943]])

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

>>> qmc.discrepancy(sample)
0.0196... # random

样本可以缩放到边界。

>>> l_bounds = [0, 2]
>>> u_bounds = [10, 5]
>>> qmc.scale(sample, l_bounds, u_bounds)
array([[1.54532796, 3.609945 ], # random
        [8.40526909, 2.1942472 ],
        [5.2177809 , 4.80031164],
        [6.80338249, 3.08795949],
        [2.65448791, 3.83491828]])

以下是其他示例,展示了构建 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。

真实世界的示例

[9] 中,使用拉丁超立方抽样 (LHS) 策略对参数空间进行抽样,以研究流行病模型中每个参数的重要性。这种分析也称为灵敏度分析。

由于问题的维度很高 (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 倍。