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()

上記のコードは、値がfloatで、合計が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日）

オリバー、mgilson、およびDougalの入力に感謝します。私の解決策を以下に示します。

  1. オリバーは多項分布のアイデアで非常に創造的でした
  2. 簡単に言えば、（1）特定のソリューションを他のソリューションよりも多く出力する可能性が非常に高いです。Dougalは、大数の法則の簡単なテスト/反例によって、多項解の空間分布が均一または正常ではないことを示しました。Dougalはまた、numpyの多項関数を使用することを提案しました。これにより、多くのトラブル、痛み、および頭痛を軽減できます。
  3. （2）の出力の問題を克服するために、RunFloat（）を使用して、より均一な分布であるように見えるもの（これはテストしていないので、表面的な外観にすぎません）を与えます。これは（1）と比べてどのくらいの違いがありますか？オフハンドはよくわかりません。でも、私の使用には十分です。
  4. numpyを使用しない代替方法を提供してくれたmgilsonに再度感謝します。

この編集のために私が作成したコードは次のとおりです。

編集＃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個の値（インデックスごとに1つ、整数も）は0からnの間でなければなりません。次に、ここのレシピに従ってください。

（または、 @ Dougalが役立つようにnumpy.random.multinomialを使用します）。

于 2012-07-08T04:00:34.113 に答える
2

@Oliverの多項アプローチと@mgilsonのコードの両方をそれぞれ100万回実行し、長さ3のベクトルの合計を10にして、それぞれの可能な結果が得られた回数を調べました。どちらも非常に不均一です。

（これから、インデックス作成のアプローチを示します。）

これは重要ですか？「通常は毎回異なるこのプロパティを持つ任意のベクトル」が必要かどうか、または各有効なベクトルが同じように発生するかどうかによって異なります。

もちろん、多項アプローチでは、3 3 4よりもはるかに可能性が高くなります0 0 10（4200倍の可能性があります）。mgilsonのバイアスは私にはあまり明白ではありませんが0 0 10、その順列ははるかに可能性が低いです（100万のうちそれぞれ約750倍）。最も一般的なのは1 4 5、その順列でした。理由はわかりませんが、確かに最も一般的で、次に1 3 6。通常、この構成では合計が高すぎることから始まります（期待値15）が、なぜそのように削減が機能するのかはわかりません。

可能なベクトル全体で均一な出力を取得する1つの方法は、拒否スキームです。K合計で長さのベクトルを取得するには、次のNようにします。

  1. との間でK均一かつ独立して整数要素を持つ長さのベクトルをサンプリングします。0N
  2. ベクトルの合計がになるまで繰り返しNます。

明らかに、これは非小型Kおよびの場合は非常に遅くなりNます。

別のアプローチは、すべての可能なベクトルに番号を割り当てることです。そのようなベクトルが(N + K - 1) choose (K - 1)あるので、その範囲内のランダムな整数を選択して、どれを使用するかを決定します。それらに番号を付ける合理的な方法の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です。

これは他のどの答えよりも遅いです：合計10の長さ3の場合、私は

  • 14.7使用してnp.random.multinomial、
  • mgilsonを使用して33.9us、
  • このアプローチで88.1私たち

考えられる結果の数が増えるにつれて、違いはさらに悪化すると思います。

誰かがどういうわけかこれらのベクトルにインデックスを付けるための気の利いた式を思いついた場合、それははるかに良いでしょう....

于 2012-07-08T06:27:07.493 に答える
1

N個の要素のパーティションのセットからK個のビンに均一にサンプリングする最も効率的な方法は、動的計画法アルゴリズムであるO（KN）を使用することです。複数選択（http://mathworld.wolfram.com/Multichoose.html）の可能性があるため、すべてを列挙するのは非常に遅くなります。棄却サンプリングやその他のモンテカルロ法も非常に遅い可能性があります。

多項分布からのサンプリングのように、人々が提案する他の方法は、一様分布からサンプルを抽出しません。

T（n、k）をn個の要素のk個のビンへの分割の数とすると、繰り返しを計算できます。

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 に答える
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 に答える
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

別のアプローチを提供するために、を実装し、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])
于 2012-07-08T11:50:18.870 に答える
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 に答える