Python的拉丁超立方采样

10

我希望对一个在多个维度(2、3、4)中定义的函数所表示的分布进行取样:

f(x, y, ...) = ...

分布可能很丑陋,非标准的(如3D数据样条曲线,高斯总和等)。 为此,我想在2..4维空间中均匀采样,然后用另一个随机数接受或拒绝给定点到我的样本中。

  1. 是否有现成的Python库可以实现这个目的?

  2. 是否有Python库可以使用拉丁超立方采样或其他均匀采样方法生成这个2..4维空间中的点? 独立随机数的暴力采样通常会导致空间中更密集和不那么密集的区域。

  3. 如果1)和2)不存在,是否有人好心分享他的实现相同或类似问题的代码。

我将在Python代码中使用它,但也可以参考其他解决方案。

5个回答

11

我想这个回答有些晚,但也是为了未来的访问者。我刚刚在git上发布了 多维均匀拉丁超立方抽样的实现 。它很简单,易于使用。如果您的变量是独立的,您可以使用拉丁超立方抽样生成在n维中采样的均匀随机变量。下面是一个示例图,比较了具有零相关性的两个维度中的蒙特卡罗和多维均匀拉丁超立方(LHS-MDU)抽样。

import lhsmdu
import matplotlib.pyplot as plt
import numpy

l = lhsmdu.sample(2,10) # Latin Hypercube Sampling of two variables, and 10 samples each.
k = lhsmdu.createRandomStandardUniformMatrix(2,10) # Monte Carlo Sampling

fig = plt.figure()
ax = fig.gca()
ax.set_xticks(numpy.arange(0,1,0.1))
ax.set_yticks(numpy.arange(0,1,0.1))
plt.scatter(k[0], k[1], color="b", label="LHS-MDU")
plt.scatter(l[0], l[1], color="r", label="MC")
plt.grid()
plt.show()

MCS versus LHS


这对我来说立即失败了:ndexError: only integers, slices (:), ellipsis (...), numpy.newaxis (None) and integer or boolean arrays are valid indices. 失败的行是 l = lhsmdu.sample(2, 10)(我在使用Python 3)。 - Santi Peñate-Vera
@SantiPeñate-Vera 很抱歉回复晚了。目前这个只支持Python 2.7,如果你还打算使用它,我可以花些时间使其兼容Python 3。 - Sahil M
我现在没有使用该库的紧迫性,但如果您将其移植到Python 3,我一定会使用它。 - Santi Peñate-Vera
好的。假设你已经正确安装了所有依赖项,它应该可以工作了。你能检查一下吗? - Sahil M
1
我尝试使用pip安装“lhsmdu”,但是没有找到这个名称的模块。 - Alexander Caskie
2
之前它不在pip上。我现在已经上传了。你能再检查一遍吗? - Sahil M

7

从1.7版本开始,拉丁超立方抽样已成为SciPy的一部分。请参见文档

from scipy.stats.qmc import LatinHypercube

engine = LatinHypercube(d=2)
sample = engine.random(n=100)

它支持居中、强度和优化。


6

现在,pyDOE库提供了一种生成基于拉丁超立方体的样本的工具。

https://pythonhosted.org/pyDOE/randomized.html

要在n个维度上生成样本:

lhs(n, [samples, criterion, iterations])

其中n是维度的数量,样本是样本空间的总数。


2
我该如何指定种子值,以便每次都能获得相同的结果? - Realhermit
不确定是否相关,但我将这些值保存在文本文件中,以便以后使用。就个人而言,我没有检查过结果是否在尝试之间有所变化,但认为在硬盘驱动器上备份一份副本可能是个好主意。 - Qiangzini
问题在于我事先不知道我的“n”和“p”参数。然而,我能够在不修改pyDOE源代码的情况下解决我的问题,方法是在导入pyDOE之后在我的主脚本中导入numpy,然后执行numpy.random.seed(seed_value)。这确保了运行之间的数字相同。 - Realhermit
pyDOE2是pyDOE的一个分支,目前仍在维护中。它的一个改进是你可以直接将随机种子传递给lhs函数。 - Nathan

3

这是 Sahil M 对 Python 3 的回答的更新(从 Python 2 升级到 Python 3 并进行了一些代码更改以匹配代码和图像):

import lhsmdu
import matplotlib.pyplot as plt
import numpy

l = lhsmdu.sample(2,10) # Latin Hypercube Sampling of two variables, and 10 samples each.
k = lhsmdu.createRandomStandardUniformMatrix(2,10) # Monte Carlo Sampling

fig = plt.figure()
ax = fig.gca()
ax.set_xticks(numpy.arange(0,1,0.1))
ax.set_yticks(numpy.arange(0,1,0.1))
plt.scatter([k[0]], [k[1]], color="r", label="MC")
plt.scatter([l[0]], [l[1]], color="b", label="LHS-MDU")
plt.legend()
plt.grid()
plt.show()

代码可视化输出

我在运行这个脚本时遇到了Python内存错误。有什么建议为什么会出现这种情况,或者如何更改脚本以避免未来再次出现此问题?


需要更多信息才能回答这个问题。您可以将其作为一个新问题或在Github上发布一个问题吗? - Sahil M
1
好的,我会做到。对于小样本,它表现得非常出色。但是对于更大的样本和/或更多的变量,有时会崩溃或生成所需的时间非常漫长(这在复杂操作方面是可以接受的,但更高的性能将更有用)。 - Matthi9000
如果有人在这个答案里看到,该算法已经更新,现在运行速度快了约20倍。 - Sahil M
感谢更新。但是,生成9D 1000个样本的速度太慢了。 - Zeel B Patel

1
这个二维示例在两个维度上均匀采样,以恒定概率选择每个点(因此保持二项分布的点数),从样本空间中随机且不重复地选择这些点,并生成一对向量,然后您可以将其传递给您的函数f:
import numpy as np
import random
resolution = 10
keepprob = 0.5
min1, max1 = 0., 1.
min2, max2 = 3., 11. 
keepnumber = np.random.binomial(resolution * resolution, keepprob,1)
array1,array2  = np.meshgrid(np.linspace(min1,max1,resolution),np.linspace(min2,max2,resolution))
randominixes  = random.sample(list(range(resolution * resolution)), int(keepnumber))
randominixes.sort()
vec1Sampled,vec2Sampled  = array1.flatten()[randominixes],array2.flatten()[randominixes]

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接