无偏返回n个随机正数（>=0）使其和==总和

Question

我想找一个算法或者建议，来改进我的代码，让它生成一组随机数，这些随机数的总和等于一个指定的数字。因为我下面的代码总是偏向于生成较大的数字，所以前面的数字通常会比较高。

有没有什么方法可以让选择数字的过程更高效呢？

#!/usr/bin/python
'''
  Generate a list of 'numbs' positive random numbers whose sum = 'limit_sum'
'''

import random


def gen_list(numbs, limit_sum):
  my_sum = []
  for index in range(0, numbs):
    if index == numbs - 1:
      my_sum.append(limit_sum - sum(my_sum))
    else:
      my_sum.append(random.uniform(0, limit_sum - sum(my_sum)))

  return my_sum

#test
import pprint
pprint.pprint(gen_list(5, 20))
pprint.pprint(gen_list(10, 200))
pprint.pprint(gen_list(0, 30))
pprint.pprint(gen_list(1, 10))

输出结果

## output

[0.10845093828525609,
 16.324799712999706,
 0.08200162072303821,
 3.4534885160590041,
 0.031259211932997744]

[133.19609626532952,
 47.464880208741029,
 8.556082341110228,
 5.7817325913462323,
 4.6342577008233716,
 0.22532341156764768,
 0.0027495225618908918,
 0.064738336208217895,
 0.028888697891734455,
 0.045250924420116689]

[]

[10]

Answer 1

这是我会怎么做的：

1. 生成 n-1 个随机数，范围在 [0,max] 之间。
2. 把这些数字进行排序。
3. 对于排序后的列表中的每一对数字，取第 i 个和第 (i+1) 个数字，创建一个区间 (i,i+1)，并计算这个区间的长度。最后一个区间从最后一个数字开始，到 max 结束，第一个区间从 0 开始，到列表中的第一个数字结束。

现在，这些区间的长度加起来总是等于 max，因为它们只是表示了 [0,max] 之间的不同部分。

代码（用 Python 写的）：

#! /usr/bin/env python
import random

def random_numbers(n,sum_to):
    values=[0]+[random.randint(0,sum_to) for i in xrange(n-1)]+[sum_to]
    values.sort()
    intervals=[values[i+1]-values[i] for i in xrange(len(values)-1)]
    return intervals

if __name__=='__main__':
    print random_numbers(5,100)

Answer 2

为什么不直接生成正确数量的均匀分布随机数，把它们加起来再调整一下呢？

补充说明一下：你想要N个数字，它们的总和是S？那么可以生成N个在[0,1)区间内均匀分布的随机数，或者你用的随机数生成器能产生的其他范围的数字。把这些数字加起来，它们的总和会是s（假设），而你想要的是总和为S，所以你可以把每个数字都乘以S/s。这样一来，这些数字就均匀随机分布在[0,S/s)这个范围内了，我想是这样的。

Answer 3

好的，我们来解决这个问题，假设我们的需求是生成一个长度为 N 的随机向量，这个向量在允许的范围内是均匀分布的，换句话说就是：

给定：

• 一个期望的长度 L，
• 一个期望的总和 S，
• 每个数值的允许范围 [0,B]，

生成一个长度为 N 的随机向量 V，使得这个随机变量 V 在它允许的空间内是均匀分布的。

---

我们可以简化这个问题，注意到我们可以计算 V = U * S，其中 U 是一个类似的随机向量，它的总和为 1，允许的数值范围是 [0,b]，其中 b = B/S。这个 b 的值必须在 1/N 和 1 之间。

---

首先考虑 N = 3。允许的值 {U} 的空间是一个平面的一部分，这个平面垂直于向量 [1 1 1]，并且通过点 [1/3 1/3 1/3]，这个平面位于一个每个分量在 0 和 b 之间的立方体内。这个点的集合 {U} 形状像一个六边形。

(待定：图片。我现在无法生成，需要访问 MATLAB 或其他可以做 3D 绘图的程序。我的 Octave 安装无法做到这一点。)

最好使用一个正交归一的加权矩阵 W（可以参考我其他的回答），其中一个向量是 [1 1 1]/sqrt(3)。这样的矩阵是：

octave-3.2.3:1> A=1/sqrt(3)
   A =  0.57735
octave-3.2.3:2> K=1/sqrt(3)/(sqrt(3)-1)
   K =  0.78868
octave-3.2.3:3> W = [A A A; A 1-K -K; A -K 1-K]
   W =

     0.57735   0.57735   0.57735
     0.57735   0.21132  -0.78868
     0.57735  -0.78868   0.21132

这个矩阵也是正交归一的（W*W = I）。

如果你考虑立方体的点 [0 0 b]、[0 b b]、[0 b 0]、[b b 0]、[b 0 0] 和 [b 0 b]，这些点形成一个六边形，并且都距离立方体的对角线 b*sqrt(2/3) 的距离。这些点不满足我们的问题，但稍后会用到。其他两个点 [0 0 0] 和 [b b b] 在立方体的对角线上。

正交归一的加权矩阵 W 使我们能够生成在 {U} 中均匀分布的点，因为正交归一矩阵是坐标变换，它们只会旋转/反射，而不会缩放或扭曲。

我们将生成在 W 定义的坐标系统中均匀分布的点。第一个分量是立方体对角线的轴。U 的分量总和完全依赖于这个轴，而与其他轴无关。因此，这个轴上的坐标被强制为 1/sqrt(3)，对应于点 [1/3, 1/3, 1/3]。

其他两个分量在与立方体对角线垂直的方向上。由于离对角线的最大距离是 b*sqrt(2/3)，我们将生成均匀分布的数字 (u,v)，范围在 -b*sqrt(2/3) 到 +b*sqrt(2/3) 之间。

这给我们一个随机变量 U' = [1/sqrt(3) u v]。然后我们计算 U = U' * W。结果中的一些点可能会超出允许的范围（U 的每个分量必须在 0 和 b 之间），如果是这样，我们就拒绝这个值并重新开始。

换句话说：

1. 生成独立的随机变量 u 和 v，它们在 -b*sqrt(2/3) 和 +b*sqrt(2/3) 之间均匀分布。
2. 计算向量 U' = [1/sqrt(3) u v]
3. 计算 U = U' * W。
4. 如果 U 的任何分量超出范围 [0,b]，就拒绝这个值并回到第 1 步。
5. 计算 V = U * S。

---

对于更高维度的情况，解决方案类似（在超立方体主对角线垂直的超平面部分内均匀分布的点）：

预先计算一个秩为 N 的加权矩阵 W。

1. 生成独立的随机变量 u₁、u₂、... u_N-1，每个在 -b*k(N) 和 +b*k(N) 之间均匀分布。
2. 计算向量 U' = [1/N u₁, u₂, ... u_N-1]
3. 计算 U = U' * W。（实际上有一些捷径，不必真的构造和乘以 W。）
4. 如果 U 的任何分量超出范围 [0,b]，就拒绝这个值并回到第 1 步。
5. 计算 V = U * S。

范围 k(N) 是一个与 N 相关的函数，表示边长为 1 的超立方体的顶点到其主对角线的最大距离。我不确定一般公式是什么，但对于 N = 3 是 sqrt(2/3)，对于 N = 5 是 sqrt(6/5)，可能在某处有公式。