在超平面上进行均匀采样

Question

给定向量大小N，我想生成一个向量<s1,s2, ..., sn>，使得s1+s2+...+sn = S。已知0<S<1和si < S。此外，所生成的向量应该是均匀分布的。任何用C编写的代码都可以帮助解释问题！

- Fshly

---

一个简单的解决方案是迭代地在范围(0, S-sum)内生成数字，其中sum是迄今为止生成的所有数字的总和，然后对列表进行洗牌。不过我觉得这个方法并不够均匀。:| - amit

一个问题：s1，s2，...sn元素是否来自输入向量？如果是的话 - 该问题等同于子集和问题，我们无法有效地知道在给定向量中是否存在一组数字总和为S，更不用说找到其中n个随机样本了。 - amit

Answer 1

这段代码在这里似乎可以解决问题，但它相当复杂。

我可能会选择一个更简单的基于拒绝的算法，即：从超平面的法向量开始，在n维空间中选择一个正交基。将每个点（S，0，0，0..0），（0，S，0，0..0）转换为该基础并沿着每个基向量存储最小值和最大值。在新的基础上均匀采样每个分量，除了第一个分量（法向量），它总是S，然后转换回原始空间并检查是否满足约束条件。如果不满足，则重新采样。

P.S. 我认为这更多是一个数学问题，实际上，询问http://maths.stackexchange.com或http://stats.stackexchange.com可能是一个好主意。

- Qnan

---

n个实数的和恰好为S的概率是0。使用浮点数时，这个概率大约是2^-50。虽然这个想法对于整数可能还可以，但对于实数来说则失败得很惨。 - amit

看起来不错，但从第一眼看起来似乎不是“拒绝采样”。您能否还添加一个简短的描述，说明如何执行此方法？那里的“抽象”并不足以理解如何实现它。 - amit

我不确定我是否理解了。您能否解释一下为什么修改后的拒绝将不会是概率0，以匹配每次迭代中的总和（s1 + s2 + s3 + ... + sn）？对于足够大的n，它是正态分布的，而正态分布数值成为任意数字的概率是0。 - amit

我不打算匹配总和，而是只采样 n-1 个维度。这样，在变换后，任何结果点都将在超平面上，但我们仍然需要检查它是否在界限内。例如，对于三维空间，我们将从超平面上的某个正方形区域进行抽样，然后拒绝不落入三角形内部的样本。 - Qnan

它的时间复杂度是多少？越高效越好。 - Fshly

@Fshly 每个样本的生成时间复杂度为O(n)，其中n是维度。但是随着n的增长，接受概率会下降（通常在基于拒绝的方法中都会如此）。如果选择最优的基底，接受概率将等于从超平面采样到其最小包围的n-1维立方体的比例。我上面提到的Matlab代码对于较高的n值更有效，并且重新实现为C语言应该不会花费太长时间，因为它似乎没有使用比矩阵运算更复杂的东西。 - Qnan

Answer 2

问题可以映射到线性多面体的采样问题上，常见的方法有蒙特卡罗方法、随机游走和hit-and-run方法（参见https://www.jmlr.org/papers/volume19/18-158/18-158.pdf进行简短比较）。它与线性规划有关，并可扩展到流形。

此外，在组成数据分析中还有对多面体的分析，例如https://link.springer.com/content/pdf/10.1023/A:1023818214614.pdf，提供了平面和多面体之间可逆变换，可用于采样。

如果您正在处理低维度，则还可以使用拒绝采样。这意味着您首先在包含多面体的平面上进行采样（由不等式定义）。后一种方法易于实现（当然是浪费的），下面是GNU Octave（我让问题的作者重新实现为C）的示例代码。

第一个要求是获得垂直于超平面的向量。对于N个变量的总和，这是n = (1,...,1)。第二个要求是平面上的一个点。对于您的示例，可以是p = (S,...,S)/N。

现在，平面上的任何点都满足n^T * (x - p) = 0

我们还假设x_i >= 0

有了这些条件，您可以在平面上计算正交基（向量n的零空间），然后在该基础上创建随机组合。最后，您将映射回原始空间并对生成的样本应用约束。

# Example in 3D
dim = 3;
S   = 1;
n   = ones(dim, 1); # perpendicular vector
p   = S * ones(dim, 1) / dim;
 
# null-space of the perpendicular vector (transposed, i.e. row vector)
# this generates a basis in the plane
V = null (n.');

# These steps are just to reduce the amount of samples that are rejected
# we build a tight bounding box
bb = S * eye(dim); # each column is a corner of the constrained region
# project on the null-space
w_bb = V \ (bb - repmat(p, 1, dim));
wmin = min (w_bb(:));
wmax = max (w_bb(:));

# random combinations and map back
nsamples = 1e3;
w        = wmin + (wmax - wmin) * rand(dim - 1, nsamples);
x        = V * w + p;

# mask the points inside the polytope
msk = true(1, nsamples);
for i = 1:dim
  msk &= (x(i,:) >= 0);
endfor

x_in = x(:, msk); # inside the polytope (your samples)
x_out = x(:, !msk); # outside the polytope

# plot the results
scatter3 (x(1,:), x(2,:), x(3,:), 8, double(msk), 'filled');
hold on
plot3(bb(1,:), bb(2,:), bb(3,:), 'xr')
axis image

- JuanPi

Answer 3

可能的一个想法是：在某个包含体积中生成许多均匀分布的点，并将它们投射到平面的目标部分上。

为了获得均匀分布，该体积必须呈现出平面部分的形状，但沿着平面法线添加边缘。

要在这样的体积中均匀地生成点，我们可以将其包围在一个立方体中，并拒绝所有位于体积之外的点。

1. 选择边缘，让我们以简单起见取margin=S（一旦margin为正值，它只会影响性能）
2. 在立方体[-M，S+M] x [-M，S+M] x [-M，S+M]中生成一个点
3. 如果到平面的距离大于M，则拒绝该点并返回#2
4. 将点投影到平面上
5. 检查投影是否落在[0，S] x [0，S] x [0，S]内，否则 - 拒绝并返回＃2
6. 将此点添加到结果集中，如果需要更多点则返回＃2

- maxim1000