在一个范围内生成 N 个具有恒定总和的随机数

Question

我想从 [a,b] 之间的特定分布（例如均匀随机）中生成 N 个随机数，总和为常数 C。我尝试了一些我自己能想到的解决方案，其中一些提出了类似的线程，但它们中的大多数要么适用于有限形式的问题，要么我无法证明结果仍然遵循所需的分布。

我尝试过的： 生成 N 个随机数，将它们除以它们的总和，然后乘以所需的常数。这似乎可行，但结果不遵循数字应在 [a:b] 内的规则。

生成 N-1 个随机数加上 0 和所需的常数 C 并对它们进行排序。然后计算每两个连续数字之间的差异，差异就是结果。这再次与 C 相加，但与最后一个方法有相同的问题（范围可以大于 [a:b]。

我还尝试生成随机数，并始终以保持所需总和和范围的方式跟踪最小值和最大值，并提出以下代码：

bool generate(function<int(int,int)> randomGenerator,int min,int max,int len,int sum,std::vector<int> &output){
    /**
    * Not possible to produce such a sequence
    */
if(min*len > sum)
    return false;
if(max*len < sum)
    return false;

int curSum = 0;
int left = sum - curSum;
int leftIndexes = len-1;
int curMax = left - leftIndexes*min;
int curMin = left - leftIndexes*max;

for(int i=0;i<len;i++){
    int num = randomGenerator((curMin< min)?min:curMin,(curMax>max)?max:curMax);
    output.push_back(num);
    curSum += num;
    left = sum - curSum;
    leftIndexes--;
    curMax = left - leftIndexes*min;
    curMin = left - leftIndexes*max;
}

return true;
}

这似乎可行，但结果有时非常不准确，我认为它不遵循原始分布（例如统一）。例如：

//10 numbers within [1:10] which sum to 50:
generate(uniform,1,10,10,50,output);
//result:
2,7,2,5,2,10,5,8,4,5 => sum=50
//This looks reasonable for uniform, but let's change to 
//10 numbers within [1:25] which sum to 50:
generate(uniform,1,25,10,50,output);
//result:
24,12,6,2,1,1,1,1,1,1 => sum= 50

注意输出中有多少个。这听起来可能是合理的，因为范围更大。但它们看起来并不像均匀分布。 我不确定即使有可能实现我想要的，也可能是限制因素使问题无法解决。

Answer 1

虽然这是一个老话题，但我想我有一个想法。考虑我们想要 N 个随机数，总和为 C，并且每个随机数介于 a 和 b 之间。为了解决问题，我们创建了 N 个洞并准备了 C 个球，每次我们问每个洞“你想要另一个球吗？”。如果没有，我们传到下一个洞，否则，我们把球放进洞里。每个洞都有一个上限值：b-a。如果某个洞达到上限值，则始终传递到下一个洞。

示例：
0 到 2 之间的 3 个随机数，总和为 5。

模拟结果：
第一次运行：-+-
第二次运行：++-
第三次运行：---
第 4 次运行：+*+
最终：221

-:垃圾球
+:接球
*:全通

Answer 2

如果您希望样本遵循均匀分布，则问题会简化为生成总和 = 1 的 N 个随机数。这又是 Dirichlet 分布的一种特殊情况，但也可以使用指数分布。方法如下：

1. 取一个均匀的样本 v₁ … v_N，所有 v_i 都在 0 和 1 之间。
2. 对于所有 i，1i := -ln v_i（注意 u_i > 0)。
3. 将 u_i 归一化为 p_i := u_i/s 其中 s 是总和 u₁+...+u_N.

p₁..p_N 是均匀分布的（在 dim N-1 的单纯形中），它们的和为 1。

您现在可以将这些 p_i 乘以您想要的常数 C，然后通过将其他一些常数 A 相加来转换它们

q_i := A + p_i*C.

编辑 3

为了解决 cmets 中提出的一些问题，让我添加以下内容：

• 为确保最终的随机序列落在区间 [a,b] 中，选择上面的常数 A 和 C 为 A := a 和 C := ba，即取 q_i = a + p_i*(ba)。由于 p_i 在 (0,1) 范围内，所有 q_i 将在 [a,b] 范围内。
• 如果 v_i 恰好为 0，则不能取（负）对数 -ln(v_i)，因为 ln() 未定义为 0。概率这种事件的发生率极低。但是，为了确保不发出错误信号，上述第 1 项中 v₁ ... v_N 的生成必须以特殊方式威胁 0 的任何出现：将 -ln(0) 视为 +infinity（记住：当 x->0 时，ln(x) -> -infinity）。因此总和 s = +infinity，这意味着 p_i = 1 和所有其他 p_j = 0。如果没有这个约定，序列 (0...1.. .0) 永远不会生成（非常感谢@Severin Pappadeux 的这个有趣的评论。）
• 正如@Neil Slater 对问题所附的第 4 条评论中所解释的，在逻辑上不可能满足原始框架的所有要求。因此，任何解决方案都必须将约束放松到原始约束的适当子集。 @Behrooz 的其他 cmets 似乎证实这在这种情况下就足够了。

编辑 2

在 cmets 中又提出了一个问题：

为什么重新调整一个统一的样本是不够的？

换句话说，我为什么要费心去取负对数？

原因是，如果我们只是重新缩放，那么生成的样本将不会均匀分布在片段 (0,1)（或最终样本的 [a,b] 中。）

为了形象化，让我们考虑 2D，即，让我们考虑 N=2 的情况。一个均匀的样本 (v₁,v₂) 对应于正方形中的一个随机点，原点 (0,0) 和角点 (1,1)。现在，当我们将这样一个点除以总和 s=v₁+v₂ 进行归一化时，我们所做的是将该点投影到对角线上，如图所示（请记住，对角线是线 x + y = 1）：

但是考虑到更靠近从 (0,0) 到 (1,1) 的主对角线的绿色线比靠近 x 轴和 y 轴的橙色线长，因此投影往往会累积更多围绕投影线的中心（蓝色），缩放样本所在的位置。这表明简单的缩放不会在所描绘的对角线上产生均匀的样本。另一方面，可以在数学上证明负对数确实产生了所需的均匀性。因此，我不会复制粘贴数学证明，而是邀请所有人实现这两种算法并检查结果图的行为是否与此答案描述的一样。

（注意：here 是一篇关于这个有趣主题的博文，并应用于石油和天然气行业）

Answer 3

对于我的回答，我假设我们有一个均匀分布。

由于我们有一个均匀分布，C 的每个元组都有相同的发生概率。例如对于a = 2, b = 2, C = 12, N = 5，我们有15 可能的元组。其中10 以2 开头，4 以3 开头，1 以4 开头。这给出了从1 到15 中选择一个随机数以选择第一个元素的想法。从1 到10，我们选择2，从11 到14，我们选择3，对于15，我们选择4。然后我们继续递归。

#include <time.h>
#include <random>

std::default_random_engine generator(time(0));
int a = 2, b = 4, n = 5, c = 12, numbers[5];

// Calculate how many combinations of n numbers have sum c
int calc_combinations(int n, int c) {
    if (n == 1) return (c >= a) && (c <= b);
    int sum = 0;
    for (int i = a; i <= b; i++) sum += calc_combinations(n - 1, c - i);
    return sum;
}

// Chooses a random array of n elements having sum c
void choose(int n, int c, int *numbers) {
    if (n == 1) { numbers[0] = c; return; }

    int combinations = calc_combinations(n, c);
    std::uniform_int_distribution<int> distribution(0, combinations - 1);
    int s = distribution(generator);
    int sum = 0;
    for (int i = a; i <= b; i++) {
        if ((sum += calc_combinations(n - 1, c - i)) > s) {
            numbers[0] = i;
            choose(n - 1, c - i, numbers + 1);
            return;
        }
    }
}

int main() { choose(n, c, numbers); }

可能的结果：

2
2
3
2
3

由于组合计算中的溢出（除非我们使用大整数库）、此计算所需的时间以及需要任意大的随机数，因此该算法无法很好地扩展 N。

Answer 4

好吧，对于 n=10000，我们不能有一个小数字不是随机的吗？

可能会生成序列直到达到sum > C-max，然后只需输入一个简单的数字即可。

万分之一更像是系统中非常小的噪音。

Answer 5

让我们尝试简化问题。 通过减去下界，我们可以将其简化为在 [0,ba] 中找到 N 个数，使得它们的和为 C-Na。

重命名参数，我们可以在[0,m]中寻找N个数，其和为S。

现在的问题类似于将长度为 S 的片段划分为 N 个长度为 [0,m] 的不同子片段。

我认为这个问题根本无法解决。

如果 S=1、N=1000 和 m 大于 0，则唯一可能的重新分区是一个 1 和 999 个零，这与随机分布完全不同。

N、m和S之间存在相关性，即使选取随机值也不会使其消失。

对于最均匀的重新划分，子段的长度将遵循平均值为S/N的高斯曲线。

如果你以不同的方式调整你的随机数，你最终会得到任何偏差，但最终你永远不会同时拥有统一的 [a,b] 重新分区和 C 的总长度，除非你的 [a ,b] 区间恰好是 2C/Na。