用numpy生成拉丁超立方体样本
我编写了一些代码来生成拉丁超立方体样本，用于高维参数空间上的插值。拉丁超立方体本质上是超立方体上点的集合，这些点被放置在立方体/矩形网格上，它具有两个点不共享任何单个坐标的特性，并且每一行/列/高维轴都被采样一次。例如，对于2个维度和4个总样本，这是一个拉丁超立方体：
| | | |x|
| |x| | |
|x| | | |
| | |x| |产生这一结果的一个算法是生成从0到N1的整数的D随机排列，其中D是维数，N是期望的样本数。将它们合并到一个DxN矩阵中，会给出一个坐标列表，它将形成一个拉丁超立方体。
我还想生成一种特殊类型的拉丁超立方体，称为对称拉丁超立方体。这个性质实质上意味着超立方体在空间反演下是不变的。表达这个条件的一个方法是(如果样本是零索引的)如果超立方体有样本(i，j，.，k)，对于某些整数指数i，j，k，则它也有样本(n-1-i，n-1-j，…，n-1-k)，其中n是样本数。对于2D和n=4，这是对称拉丁超立方体的一个例子：
| |x| | |
| | | |x|
|x| | | |
| | |x| |我要这么做的原因是，我想要计算一个计算成本很高的函数，它依赖于许多参数。我在一组点上为这个函数预生成一个值列表，并使用参数上的插值将这个列表扩展到一个集合范围内的任何想要的点。拉丁超立方体样本用于预计算，而不是网格，其思想是更有效地对函数的行为进行采样，从而降低插值误差。
我已经测试了下面的代码，使用了自动化测试，并对普通2D案例的结果进行了可视化检查。高维案例的统计检验使我确信，它在那里也是成功的。该代码可以生成随机拉丁超立方体或非一般对称超立方体(我确信我的第二个算法原则上不能生成任何可能的对称拉丁超立方体，只是一个特定类型)。而且，对称算法对于奇数样本并不能正常工作，所以我完全阻止了这些样本的生成。
我的代码可以与规范和性能测试这里一起看到，尽管目前还没有可视化测试。
我知道我的项目结构对于这里的代码数量来说可能是过分的。我主要是寻找任何关于我的算法的洞察力，因为虽然它们非常快，即使是大尺寸和样本大小，我不怀疑存在改进。不过，欢迎任何意见。
import numpy as np
class LatinSampler:
def _check_num_is_even(self, num):
if num % 2 != 0:
raise ValueError("Number of samples must be even")
def get_lh_sample(self, param_mins, param_maxes, num_samples):
dim = param_mins.size
latin_points = np.array([np.random.permutation(num_samples) for i in range(dim)]).T
lengths = (param_maxes - param_mins)[None, :]
return lengths*(latin_points + 0.5)/num_samples + param_mins[None, :]
def get_sym_sample(self, param_mins, param_maxes, num_samples):
self._check_num_is_even(num_samples)
dim = param_mins.size
even_nums = np.arange(0, num_samples, 2)
permutations = np.array([np.random.permutation(even_nums) for i in range(dim)])
inverses = (num_samples - 1) - permutations
latin_points = np.concatenate((permutations,inverses), axis=1).T
lengths = (param_maxes - param_mins)[None, :]
return lengths*(latin_points + 0.5)/num_samples + param_mins[None, :]编辑:我忽略了这一点，但是这段代码本身并不返回拉丁超级立方体，而是拉长/压缩和移位拉丁超立方体，这样点就可以超过任意范围的值。但是，如果省略最后一个操作(由lengths和add param_mins乘以)，核心代码将生成一个拉丁超立方体。
回答 2
Code Review用户
发布于 2019-07-06 19:58:33
一旦我理解了您想要完成的任务，您的代码对我来说就变得更加清晰了，我可以提供一些想法。
关于对称版本，您实际上只是从可用的对称拉丁超立方体的子集中取样。参数空间的一半永远不会被选择作为样本。考虑使用n=4的2D示例：
even_nums = [0, 2]，因此，A，B，C，D是用这种方法可以生成的4个可能的样本。事实上，有两个可能的拉丁超级立方体：A，D和B、C。A'，B'，C'，D‘表示它们各自的逆位置，然后加起来。正如你所看到的，一半的字段是空的，永远不会被选择作为一个样本。
与置换索引不同，我将直接对对称点对进行采样，因为它允许您在实现均匀采样的同时控制约束：
import random
import numpy as np
num_samples = 5
dim = 2
available_indices = [set(range(num_samples)) for _ in range(dim)]
samples = []
# if num_samples is odd, we have to chose the midpoint as a sample
if num_samples % 2 != 0:
k = num_samples//2
samples.append([k] * dim)
for idx in available_indices:
idx.remove(k)
# sample symmetrical pairs
for _ in range(num_samples//2):
sample1 = list()
sample2 = list()
for idx in available_indices:
k = random.sample(idx, 1)[0]
sample1.append(k)
sample2.append(num_samples-1-k)
idx.remove(k)
idx.remove(num_samples-1-k)
samples.append(sample1)
samples.append(sample2)
samples = np.array(samples)Code Review用户
发布于 2021-08-01 12:30:37
您可以使用我在1.7版中添加到LatinHypercube中的SciPy类。这可以替换这里的大部分代码。
要生成LHS，只需生成一个统一的样本U(0，1)，然后用整数的随机排列将其移动，如下所示：
import numpy as np
n, d = ...
rng = np.random.default_rng()
samples = rng.uniform(size=(n, d))
perms = np.tile(np.arange(1, n + 1), (d, 1))
for i in range(d):
rng.shuffle(perms[i, :])
perms = perms.T
samples = (perms - samples) / nhttps://codereview.stackexchange.com/questions/223569
复制相似问题
腾讯云开发者
