8

创建一个总和为 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，也发布了一个很好的问题解释。

4

7 回答 7

3

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

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

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

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

（或者使用numpy.random.multinomial作为@Dougal 有用的建议）。

于 2012-07-08T04:00:34.113 回答
2

我刚刚将@Oliver的多项式方法和@mgilson 的代码分别运行了一百万次，长度为3 的向量总和为10，并查看了每个可能结果出现的次数。两者都非常不均匀：

（我将展示索引方法。）

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

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

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

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

显然，对于非小型K和N.

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

请注意，K向量的最后一个 (th) 元素由第一个元素的和唯一确定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 我们采用这种方法

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

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

于 2012-07-08T06:27:07.493 回答
1

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

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)
于 2012-07-08T03:56:50.110 回答
1

从 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）中采样一组隐藏状态。

于 2012-07-08T15:27:22.590 回答
0

只是为了给你另一种方法，实现 apartition_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])
于 2012-07-08T11:50:18.870 回答
0

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

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)]
于 2012-07-08T09:34:29.890 回答
0

用什么：

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]
于 2015-06-02T13:56:27.977 回答