如何创建总和为x的随机整数向量列表

Question

创建一个总和为 X（例如 X=1000）的随机向量非常简单：

import random
def RunFloat():
    Scalar = 1000
    VectorSize = 30
    RandomVector = [random.random() for i in range(VectorSize)]
    RandomVectorSum = sum(RandomVector)
    RandomVector = [Scalar*i/RandomVectorSum for i in RandomVector]
    return RandomVector
RunFloat()

上面的代码创建了一个向量，其值为浮点数，总和为 1000。

我很难创建一个简单的函数来创建一个向量，其值为整数且总和为 X（例如 X=1000*30）

import random
def RunInt():
    LowerBound = 600
    UpperBound = 1200
    VectorSize = 30
    RandomVector = [random.randint(LowerBound,UpperBound) for i in range(VectorSize)]
    RandomVectorSum = 1000*30
    #Sanity check that our RandomVectorSum is sensible/feasible
    if LowerBound*VectorSize <= RandomVectorSum and RandomVectorSum <= UpperBound*VectorSum:
        if sum(RandomVector) == RandomVectorSum:
            return RandomVector
        else:
            RunInt()  

有没有人有任何改进这个想法的建议？我的代码可能永远不会完成或遇到递归深度问题。

编辑（2012 年 7 月 9 日）

感谢 Oliver、mgilson 和 Dougal 的投入。我的解决方案如下所示。

1. Oliver 的多项式分布理念非常有创意
2. 简而言之，(1) 很可能输出某些解决方案而不是其他解决方案。 Dougal 通过大数定律的简单测试/反例证明了多项解空间分布不均匀或不正态。 Dougal 还建议使用 numpy 的多项式函数，它为我省去了很多麻烦、痛苦和头痛。
3. 为了克服 (2) 的输出问题，我使用 RunFloat() 将出现的内容（我没有对此进行测试，所以它只是一个表面的外观）变得更加均匀分布。与（1）相比，这有多大不同？我真的不知道副手。不过对我来说已经足够了。
4. 再次感谢 mgilson 提供不使用 numpy 的替代方法。

这是我为此编辑编写的代码：

编辑 #2（2012 年 7 月 11 日）

我意识到正态分布没有正确实现，我已将其修改为以下内容：

import random
def RandFloats(Size):
    Scalar = 1.0
    VectorSize = Size
    RandomVector = [random.random() for i in range(VectorSize)]
    RandomVectorSum = sum(RandomVector)
    RandomVector = [Scalar*i/RandomVectorSum for i in RandomVector]
    return RandomVector

from numpy.random import multinomial
import math
def RandIntVec(ListSize, ListSumValue, Distribution='Normal'):
    """
    Inputs:
    ListSize = the size of the list to return
    ListSumValue = The sum of list values
    Distribution = can be 'uniform' for uniform distribution, 'normal' for a normal distribution ~ N(0,1) with +/- 5 sigma  (default), or a list of size 'ListSize' or 'ListSize - 1' for an empirical (arbitrary) distribution. Probabilities of each of the p different outcomes. These should sum to 1 (however, the last element is always assumed to account for the remaining probability, as long as sum(pvals[:-1]) <= 1).  
    Output:
    A list of random integers of length 'ListSize' whose sum is 'ListSumValue'.
    """
    if type(Distribution) == list:
        DistributionSize = len(Distribution)
        if ListSize == DistributionSize or (ListSize-1) == DistributionSize:
            Values = multinomial(ListSumValue,Distribution,size=1)
            OutputValue = Values[0]
    elif Distribution.lower() == 'uniform': #I do not recommend this!!!! I see that it is not as random (at least on my computer) as I had hoped
        UniformDistro = [1/ListSize for i in range(ListSize)]
        Values = multinomial(ListSumValue,UniformDistro,size=1)
        OutputValue = Values[0]
    elif Distribution.lower() == 'normal':
        """
        Normal Distribution Construction....It's very flexible and hideous
        Assume a +-3 sigma range.  Warning, this may or may not be a suitable range for your implementation!
        If one wishes to explore a different range, then changes the LowSigma and HighSigma values
        """
        LowSigma    = -3#-3 sigma
        HighSigma   = 3#+3 sigma
        StepSize    = 1/(float(ListSize) - 1)
        ZValues     = [(LowSigma * (1-i*StepSize) +(i*StepSize)*HighSigma) for i in range(int(ListSize))]
        #Construction parameters for N(Mean,Variance) - Default is N(0,1)
        Mean        = 0
        Var         = 1
        #NormalDistro= [self.NormalDistributionFunction(Mean, Var, x) for x in ZValues]
        NormalDistro= list()
        for i in range(len(ZValues)):
            if i==0:
                ERFCVAL = 0.5 * math.erfc(-ZValues[i]/math.sqrt(2))
                NormalDistro.append(ERFCVAL)
            elif i ==  len(ZValues) - 1:
                ERFCVAL = NormalDistro[0]
                NormalDistro.append(ERFCVAL)
            else:
                ERFCVAL1 = 0.5 * math.erfc(-ZValues[i]/math.sqrt(2))
                ERFCVAL2 = 0.5 * math.erfc(-ZValues[i-1]/math.sqrt(2))
                ERFCVAL = ERFCVAL1 - ERFCVAL2
                NormalDistro.append(ERFCVAL)  
            #print "Normal Distribution sum = %f"%sum(NormalDistro)
            Values = multinomial(ListSumValue,NormalDistro,size=1)
            OutputValue = Values[0]
        else:
            raise ValueError ('Cannot create desired vector')
        return OutputValue
    else:
        raise ValueError ('Cannot create desired vector')
    return OutputValue
#Some Examples        
ListSize = 4
ListSumValue = 12
for i in range(100):
    print RandIntVec(ListSize, ListSumValue,Distribution=RandFloats(ListSize))

上面的代码可以在github上找到。这是我为学校建造的课程的一部分。 user1149913，也发布了一个很好的问题解释。

Answer 1

用什么：

import numpy as np
def RunInt(VectorSize, Sum):
    a = np.array([np.random.rand(VectorSize)])
    b = np.floor(a/np.sum(a)*Sum) 
    for i in range(int(Sum-np.sum(b))):
        b[0][np.random.randint(len(b[0]))] += 1
    return b[0]

Answer 2

从 N 个元素的分区集中均匀采样到 K 个 bin 中的最有效方法是使用动态规划算法，即 O(KN)。有多种选择 (http://mathworld.wolfram.com/Multichoose.html) 的可能性，因此枚举每一个会​​非常慢。拒绝抽样和其他蒙特卡罗方法也可能非常缓慢。

人们提出的其他方法（例如从多项式中抽样）不会从均匀分布中抽取样本。

令 T(n,k) 为将 n 个元素划分为 k 个 bin 的数量，然后我们可以计算递归

T(n,1)=1 \forall n>=0
T(n,k)=\sum_{m<=n} T(n-m,k-1)

要对总和为 N 的 K 个元素进行采样，请从循环中“向后”的 K 个多项式分布中采样： 编辑：以下多项式中的 T 在抽取每个样本之前应归一化为总和为 1。

n1 = multinomial([T(N,K-1),T(N-1,K-1),...,T(0,K-1)])
n2 = multinomial([T(N-n1,K-1),T(N-n1-1,K-1),...,T(0,K-1)])
...
nK = multinomial([T(N-sum([n1,...,n{k-1}]),1),T(N-sum([n1,...,n{k-1}])-1,1),...,T(0,1)])

注意：我允许对 0 进行采样。

此过程类似于从分段半马尔可夫模型 (http://www.gatsby.ucl.ac.uk/%7Echuwei/paper/icml103.pdf) 中采样一组隐藏状态。

Answer 3

这个版本会给出一个统一的分布：

from random import randint

def RunInt(VectorSize, Sum):
   x = [randint(0, Sum) for _ in range(1, VectorSize)]
   x.extend([0, Sum])
   x.sort()
   return [x[i+1] - x[i] for i in range(VectorSize)]

Answer 4

只是给你另一种方法，实现 partition_function(X) 并随机选择一个介于 0 和 partition_function(1000) 长度之间的数字，然后你就有了。现在您只需要找到一种有效的方法来计算配分函数。这些链接可能会有所帮助：

http://code.activestate.com/recipes/218332-generator-for-integer-partitions/

http://oeis.org/A000041

编辑： 这是一个简单的代码：

import itertools
import random
all_partitions = {0:set([(0,)]),1:set([(1,)])}

def partition_merge(a,b):
    c = set()
    for t in itertools.product(a,b):
        c.add(tuple(sorted(list(t[0]+t[1]))))
    return c

def my_partition(n):
    if all_partitions.has_key(n):
        return all_partitions[n]
    a = set([(n,)])
    for i in xrange(1,n/2+1):
        a = partition_merge(my_partition(i),my_partition(n-i)).union(a)
    all_partitions[n] = a
    return a

if __name__ == '__main__':
    n = 30
    # if you have a few years to wait uncomment the next line
    # n = 1000
    a = my_partition(n)
    i = random.randint(0,len(a)-1)
    print(list(a)[i])

Answer 5

我只是将@Oliver's multinomial approach 和@mgilson's code 分别运行了一百万次，长度为 3 的向量总和为 10，并查看了每个可能结果出现的次数。两者都非常不统一：

（我即将展示索引方法。）

这有关系吗？取决于您是否想要“具有此属性的任意向量，通常每次都不同”与每个有效向量的可能性相同。

在多项式方法中，3 3 4 的可能性当然会比 0 0 10 高得多（事实证明，可能性是 4200 倍）。 mgilson 的偏见对我来说不太明显，但0 0 10 及其排列是迄今为止最不可能的（每百万次只有约 750 次）；最常见的是1 4 5 及其排列；不知道为什么，但它们肯定是最常见的，其次是1 3 6。在这个配置中，它通常会从一个太高的总和开始（预期 15），但我不确定为什么会这样减少......

在可能的向量上获得统一输出的一种方法是拒绝方案。要获得长度为 K 且总和为 N 的向量，您需要：

1. 在0 和N 之间均匀且独立地对长度为K 的向量进行采样。
2. 重复直到向量的总和为N。

显然，对于非小型 K 和 N，这将非常缓慢。

另一种方法是为所有可能的向量分配一个编号；有(N + K - 1) choose (K - 1) 这样的向量，所以只需在该范围内选择一个随机整数来决定你想要哪个。对它们进行编号的一种合理方法是按字典顺序排列：(0, 0, 10), (0, 1, 9), (0, 2, 8), (0, 3, 7), ...。

请注意，向量的最后一个 (Kth) 元素由第一个 K-1 之和唯一确定。

我确信有一个很好的方法可以立即跳转到此列表中的任何索引，但我现在想不出来......列举可能的结果并遍历它们会起作用，但可能会比必要的慢。这是一些代码（尽管我们实际上在这里使用了反向字典排序......）。

from itertools import islice, combinations_with_replacement
from functools import reduce
from math import factorial
from operator import mul
import random

def _enum_cands(total, length):
    # get all possible ways of choosing 10 of our indices
    # for example, the first one might be  0000000000
    # meaning we picked index 0 ten times, for [10, 0, 0]
    for t in combinations_with_replacement(range(length), 10):
        cand = [0] * length
        for i in t:
            cand[i] += 1
        yield tuple(cand)

def int_vec_with_sum(total, length):
    num_outcomes = reduce(mul, range(total + 1, total + length)) // factorial(length - 1)
    # that's integer division, even though SO thinks it's a comment :)
    idx = random.choice(range(num_outcomes))
    return next(islice(_enum_cands(total, length), idx, None))

如上面的直方图所示，这实际上对于可能的结果是一致的。它也很容易适应任何单个元素的上限/下限；只需将条件添加到_enum_cands。

这比其他任何一个答案都慢：对于 sum 10 length 3，我得到

• 14.7 我们使用np.random.multinomial，
• 33.9 我们使用 mgilson 的，
• 88.1 我们采用这种方法

我预计随着可能结果数量的增加，差异会变得更糟。

如果有人想出一个漂亮的公式来索引这些向量，那就更好了......

Answer 6

我建议不要递归地这样做：

当您递归采样时，第一个索引中的值具有更大的可能范围，而后续索引中的值将受第一个值的约束。这将产生类似于exponential distribution 的东西。

相反，我建议从multinomial distribution 中取样。这将平等对待每个索引，约束总和，强制所有值为整数，并从遵循这些规则的所有可能配置中统一采样（注意：可能以多种方式发生的配置将通过它们可能发生的方式数加权）。

为了帮助您将问题与多项式表示法合并，总和为 n（整数），因此每个 k 值（每个索引一个，也是整数）必须介于 0 和 n 之间。然后按照食谱here。

（或者按照@Dougal 的建议使用numpy.random.multinomial）。

Answer 7

这是一个非常简单的实现。

import random
import math

def randvec(vecsum, N, maxval, minval):
    if N*minval > vecsum or N*maxval < vecsum:
        raise ValueError ('Cannot create desired vector')

    indices = list(range(N))
    vec = [random.randint(minval,maxval) for i in indices]
    diff = sum(vec) - vecsum # we were off by this amount.

    #Iterate through, incrementing/decrementing a random index 
    #by 1 for each value we were off.
    while diff != 0:  
        addthis = 1 if diff > 0 else -1 # +/- 1 depending on if we were above or below target.
        diff -= addthis

        ### IMPLEMENTATION 1 ###
        idx = random.choice(indices) # Pick a random index to modify, check if it's OK to modify
        while not (minval < (vec[idx] - addthis) < maxval):  #operator chaining.  If you don't know it, look it up.  It's pretty cool.
            idx = random.choice(indices) #Not OK to modify.  Pick another.

        vec[idx] -= addthis #Update that index.

        ### IMPLEMENTATION 2 ###
        # random.shuffle(indices)
        # for idx in indices:
        #    if minval < (vec[idx] - addthis) < maxval:
        #        vec[idx]-=addthis
        #        break
        #
        # in situations where (based on choices of N, minval, maxval and vecsum)
        # many of the values in vec MUST BE minval or maxval, Implementation 2
        # may be superior.

    return vec

a = randvec(1000,20,100,1)
print sum(a)