非偏差返回n个随机正数(> = 0)的列表,以便它们的总和== total_sum

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

为什么不生成正确数量的均匀分布的随机数,将它们组合起来并进行扩展？

编辑:要更清楚一点:你想要N个数加到S？因此,在区间[0,1)或RNG产生的任何内容上生成N个均匀分布的随机数.添加它们,它们将总计s(比如说),而你希望它们总计为S,所以将每个数字乘以S/s.现在我认为数字在[0,S/s]上均匀随机分布.

Answer 2

我是这样做的:

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 3

如果您正在寻找具有尽可能少的相关性的正态分布数,并且需要严格*关于此,我建议您采用以下数学方法并转换为代码.

(*严谨:其他方法的问题在于你可以在你的发行版中得到"长尾巴" - 换句话说,它很少但可能有与你预期的输出非常不同的异常值)

• 生成N-1个独立且相同分布(IID)的高斯随机变量v ₀,v ₁,v ₂,... v _N-1以匹配问题的N-1个自由度.
• 创建列向量V,其中V = [0 v ₀,v ₁,v ₂,... v _N-1 ] ^T
• 使用固定加权矩阵W,其中W由正交矩阵**组成,其顶行为[1 1 1 1 1 1 1 1 ... 1]/sqrt(N).
• 您的输出向量是乘积WV + SU/N,其中S是所需的和,U是1的列向量.换句话说,第i个输出变量=(矩阵W的行#i)和列向量V的点积加到S/N.

每个输出变量的标准偏差(我相信,现在无法验证)sqrt(N/N-1)*输入随机变量的标准差.

**正交矩阵:这是困难的部分,我在math.stackexchange.com上提出了一个问题,并且有一个简单的矩阵W可以工作,并且可以通过算法定义只有3个不同的值,所以你实际上没有构造矩阵.

W是Vw的Householder反射,其中v = [sqrt(N),0,0,0,...]和w = [1 1 1 1 1 ... 1]可以通过以下定义:

W(1,i) = W(i,1) = 1/sqrt(N)
W(i,i) = 1 - K   for i >= 2 
W(i,j) = -K      for i,j >= 2, i != j
K = 1/sqrt(N)/(sqrt(N)-1)

马克方法的问题:

为什么不生成正确数量的均匀分布的随机数,将它们组合起来并进行扩展？

如果你这样做,你得到一个"长尾"分布.这是MATLAB中的一个例子:

 >> X = rand(100000,10);
 >> Y = X ./ repmat(sum(X,2),1,10);
 >> plot(sort(Y))

我在矩阵X中生成了100,000组N = 10个数,并创建了矩阵Y,其中Y的每一行是X的相应行除以其总和(因此Y的每一行总和为1.0)

绘制Y的排序值(每个列分别排序)产生大致相同的累积分布:

真正的均匀分布将产生从0到最大值的直线.你会注意到它与真正的均匀分布有点类似,除了在有长尾的末端.在0.2和0.5之间产生了过多的数字.对于较大的N值,尾部变得更糟,因为尽管数字的平均值下降(平均值= 1/N),但最大值保持为1.0:由9个值0.0和1值1.0组成的向量是有效的并且可以通过这种方式生成,但在病理上是罕见的.

如果您不关心这一点,请继续使用此方法.并且可能有方法生成具有所需总和的"几乎" - 均匀或"几乎" - 高斯分布,这比我上面描述的更简单和更有效.但我提醒您注意并理解您选择的算法的后果.

一个没有长尾分布均匀分布的修正如下:

1. 生成向量V = N均匀分布的从0.0到1.0的随机数.
2. 找出它们的和S及其最大值M.
3. 如果S <k*M(最大值超出异常值),请返回步骤1.我不确定k用什么值,也许k = N/2？
4. 输出矢量V*S _desired/S.

MATLAB中N = 10的示例:

 >> X = rand(100000,10);
 >> Y = X ./ repmat(sum(X,2),1,10);
 >> i = sum(X,2)>(10/2)*max(X,[],2);
 >> plot(sort(Y(i,:)))

Answer 4

好吧，假设要求是生成一个长度为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}的形状像六边形。

（TBD：图片。我现在无法生成一个图像，我需要访问MATLAB或另一个可以进行3D绘图的程序。我无法安装Octave。）

最好使用一个向量= [1 1 1] / sqrt（3）的正交加权矩阵W（请参阅我的其他答案）。一种这样的矩阵是

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 bb]，[0 b 0]，[bb 0]，[b 0 0]和[b 0 b]的点，它们形成一个六边形并且都是a b * sqrt（2/3）与立方体对角线的距离。这些不能满足所讨论的问题，但是在一分钟内很有用。另外两个点[0 0 0]和[bbb]在立方体的对角线上。

正交加权矩阵W允许我们生成在{U}内均匀分布的点，因为正交矩阵是旋转/反射并且不缩放或不倾斜的坐标变换。

我们将生成在W的3个向量定义的坐标系中均匀分布的点。第一个分量是立方体对角线的轴。U分量的总和完全取决于该轴，而不完全取决于其他轴。因此，沿该轴的坐标被强制为1 / sqrt（3），它对应于点[1 / 3、1 / 3、1 / 3]。

其他两个分量的方向垂直于立方体的对角线。由于距对角线的最大距离为b * sqrt（2/3），我们将在-b * sqrt（2/3）和+ b * sqrt（2/3）之间生成均匀分布的数（u，v）。

这给了我们一个随机变量U'= [1 / sqrt（3）uv]。然后，我们计算U = U'*W。一些结果点将超出允许范围（U的每个分量必须在0到b之间），在这种情况下，我们将拒绝并重新开始。

换一种说法：

1. 生成独立的随机变量u和v，它们分别均匀地分布在-b * sqrt（2/3）和+ b * sqrt（3）之间。
2. 计算向量U'= [1 / sqrt（3）uv]
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的函数，N表示侧面1的超立方体的顶点与其主对角线之间的最大距离。我不确定通用公式，但对于N = 3是sqrt（2/3），对于N = 5是sqrt（6/5），可能在某个地方有一个公式。