scipy.stats.qmc.

拉丁超立方#

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

拉丁超立方抽样 (LHS)。

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

参数:

dint

参数空间的维度。

scramblebool，可选

如果为 False，则将样本置于多维网格的单元格内中心位置。否则，样本随机放置在网格的单元格内。

注意

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

默认值为 True。

在 1.10.0 版本中添加。

optimization{None, “random-cd”, “lloyd”}，可选

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

random-cd：坐标的随机排列，以降低中心差异。基于中心差异的最佳样本不断更新。与使用其他差异度量相比，基于中心差异的采样在 2D 和 3D 子投影方面表现出更好的空间填充稳健性。
lloyd：使用修改后的 Lloyd-Max 算法扰动样本。该过程会收敛到等距样本。

在 1.8.0 版本中添加。

在 1.10.0 版本中更改：添加 lloyd。

strength{1, 2}，可选

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

在 1.8.0 版本中添加。

rngnumpy.random.Generator，可选

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

在 1.15.0 版本中更改：作为 SPEC-007 从使用 numpy.random.RandomState 过渡到使用 numpy.random.Generator 的一部分，此关键字从 seed 更改为 rng。在过渡期间，这两个关键字将继续工作，但一次只能指定一个。在过渡期之后，使用 seed 关键字的函数调用将发出警告。经过一段弃用期后，将删除 seed 关键字。

另请参阅

准蒙特卡洛

注释

当使用 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，第 2 期：143-151，1987。

[3]

A. B. Owen，“扰乱网格求积的蒙特卡洛方差”。SIAM 数值分析杂志 34，第 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，“用于计算机实验、积分和可视化的正交数组。”Statistica Sinica，1992 年。

[8]

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

[9]

西霍尔姆，苏珊·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 倍。

方法

`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 个。
`重置`()	将引擎重置为基本状态。