以条件生成随机数列表 - numpy

Question

我希望生成一个包含15个数字的列表，这些数字之和为12，最小值为0，最大值为6。

我尝试了以下代码：

def generate(low,high,total,entity):
   while sum(entity)!=total:
       entity=np.random.randint(low, high, size=15)
   return entity

但是上述函数无法正常工作。它花费的时间太长了。 请告诉我生成这样的数字的高效方法？

- Danish

---

严格地说，它是工作的，但生成和测试通常不太高效。在这里通常需要数万次生成，才能生成一个正确的序列。 - Willem Van Onsem

@WillemVanOnsem 有没有其他更快的方法？ - Danish

Answer 1

以上内容在严格意义上是可以运行的。但对于0到6之间的15个数字，生成12的概率并不那么高。实际上，我们可以通过以下公式来计算可能性的数量：

当 0≤s≤6 时，F(s, 1) = 1 并且 F(s, n) = Σ⁶_i=0F(s-i, n-1)。

我们可以使用一个值来计算它：

from functools import lru_cache

@lru_cache()
def f(s, n, mn, mx):
    if n < 1:
        return 0
    if n == 1:
        return int(mn <= s <= mx)
    else:
        if s < mn:
            return 0
        return sum(f(s-i, n-1, mn, mx) for i in range(mn, mx+1))

这意味着在总计 4'747'561'509'943 种可能性中，有 9'483'280 种可能性生成总和为 12，占比 0.00019975%。因此，需要大约 500'624 次迭代才能得出这样的解决方案。

因此，我们最好找到一种简单直接的方法来生成这样的序列。我们可以通过每次计算生成数字的概率来实现：作为前 n 个数字之和为 s 的序列的第一个数字生成 i 的概率是 F(s-i, n-1, 0, 6)/F(s, n, 0, 6)。如果我们每次抽取一个均匀分布的数字，则这将保证我们在所有符合给定条件的值列表上生成均匀的列表，而不是与整个列表的均匀分布相匹配。

我们可以使用递归来实现此操作：

from numpy import choice

def sumseq(n, s, mn, mx):
    if n > 1:
        den = f(s, n, mn, mx)
        val, = choice(
            range(mn, mx+1),
            1,
            p=[f(s-i, n-1, mn, mx)/den for i in range(mn, mx+1)]
        )
        yield val
        yield from sumseq(n-1, s-val, mn, mx)
    elif n > 0:
        yield s

通过上述函数，我们可以生成NumPy数组：

>>> np.array(list(sumseq(15, 12, 0, 6)))
array([0, 0, 0, 0, 0, 4, 0, 3, 0, 1, 0, 0, 1, 2, 1])
>>> np.array(list(sumseq(15, 12, 0, 6)))
array([0, 0, 1, 0, 0, 1, 4, 1, 0, 0, 2, 1, 0, 0, 2])
>>> np.array(list(sumseq(15, 12, 0, 6)))
array([0, 1, 0, 0, 2, 0, 3, 1, 3, 0, 1, 0, 0, 0, 1])
>>> np.array(list(sumseq(15, 12, 0, 6)))
array([5, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1])
>>> np.array(list(sumseq(15, 12, 0, 6)))
array([0, 0, 0, 0, 4, 2, 3, 0, 0, 0, 0, 0, 3, 0, 0])

- Willem Van Onsem

---

这不是在某种奇怪的方式下定义的多项式分布吗？ - Severin Pappadeux

@SeverinPappadeux：可能是这种情况，我会尝试去看一下。我有一个想法，即(a)实现应该很简单，所以不需要“拒绝”，(2)好像你从均匀采样中抽取，然后进行拒绝。如果我看一下公式，可能是这种情况，但证明它可能需要一些工作 :) - Willem Van Onsem

如果你发现了有趣的东西，我会很感激如果你给我发送一条消息，谢谢。 - Severin Pappadeux

我更新了我的答案，认为我们提出了不同的解决方案，请看一下。 - Severin Pappadeux

Answer 2

您可以尝试以略微不同的方式实现它。

import random
def generate(low,high,goal_sum,size=15):
    output = []
    for i in range(size):
        new_int = random.randint(low,high)
        if sum(output) + new_int <= goal_sum:
            output.append(new_int)
        else:
            output.append(0)
    random.shuffle(output)
    return output

此外，如果您使用np.random.randint，您的高值实际上会是高-1。

- J Lee

Answer 3

好的，有一个简单自然的解决方案 - 使用分布，它根据定义为您提供具有固定总和的值数组。最简单的是多项式分布。唯一需要添加的代码是检查并拒绝（并重复采样），如果某个采样值超过最大值。

沿着这条线

import numpy as np

def sample_sum_interval(n, p, maxv):
    while True:
        q = np.random.multinomial(n, p)
        v = np.where(q > maxv)
        if len(v[0]) == 0: # if len(v) > 0, some values are outside the range, reject
            return q
    return None

np.random.seed(32345)

k    = 15
n    = 12
maxv = 6
p = np.full((k), np.float64(1.0)/np.float64(k), dtype=np.float64) # probabilities

q = sample_sum_interval(n, p, maxv)
print(q)
print(np.sum(q))

q = sample_sum_interval(n, p, maxv)
print(q)
print(np.sum(q))

q = sample_sum_interval(n, p, maxv)
print(q)
print(np.sum(q))

更新

我快速地查看了@WillemVanOnsem提出的方法，我认为它与我使用的多项式不同。

如果我们看一下多项式PMF，并假设所有k个数字的概率相等， p₁ = ... = p_k = 1/k，则可以将PMF写成

PMF(x₁,...x_k)=n!/(x₁!...x_k!) p₁^x₁...p_k^x_k = n!/(x₁!...x_k!) k^-x₁...k^-x_k = n!/(x₁!...x_k!) k^-Sum_ix_i = n!/(x₁!...x_k!) k^-n

显然，由于分母中的阶乘不同（当然要考虑排列），特定x1...xk组合的概率会有所不同，这与@WillemVanOnsem的方法不同，后者认为它们都有相等的出现概率，我相信。

故事的寓意是-这些方法产生不同的分布。

- Severin Pappadeux