c - 在超平面上均匀采样

Question

给定向量大小 N ，我想生成一个向量<s1,s2, ..., sn>。s1+s2+...+sn = S

已知0<S<1和si < S。生成的这些向量也应该是均匀分布的。

任何有助于解释的 C 代码都会很棒！

Answer 1

这里的代码似乎可以解决问题，尽管它相当复杂。

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

PS 我认为这更像是一个数学问题，实际上，在http://maths.stackexchange.com或http://stats.stackexchange.com上提问可能是个好主意

Answer 2

[为简单起见，我将跳过“hyper-”前缀]

一种可能的想法：在某个封闭体积中生成许多均匀分布的点，并将它们投影到平面的目标部分。

为了获得均匀分布，体积必须像平面的一部分，但沿平面法线添加边距。

为了在这样的体积中均匀地生成点，我们可以将它封闭在一个立方体中并拒绝体积之外的所有内容。

1. 选择保证金，为简单起见，我们取保证金=S（一旦保证金为正，它只会影响性能）
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 是否需要更多点

Answer 3

该问题可以映射到线性多面体上的采样问题，常见的方法是蒙特卡罗方法、随机游走和命中运行方法（参见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