生成带有条件的随机数列表 - numpy [重复]

Posted 2023-03-11

【中文标题】生成带有条件的随机数列表 - numpy [重复]

【英文标题】：Generate list of random number with condition - numpy [duplicate]

【发布时间】：2020-01-22 16:09:06

【问题描述】：

我想生成一个包含 15 个整数的列表，总和为 12，最小值为 0，最大值为 6。

我尝试了以下代码

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

但上述功能无法正常工作。这太费时间了。 请告诉我生成此类数字的有效方法？

【问题讨论】：

严格来说它可以工作，但是生成和测试通常不是很有效。在您生成正确的序列之前，它通常需要一万次生成。

@WillemVanOnsem 还有其他快捷的方法吗？

【参考方案1】：

严格来说，以上将起作用。但是对于 0 到 6 之间的 15 个数字，生成 12 个的几率并没有那么高。事实上，我们可以通过以下方式计算可能性的数量：

F(s, 1) = 1 对于 0≤s≤6 和

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

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

因此，我们应该更好地寻找一种直接的方式来生成这样的序列。我们可以通过每次计算生成一个数字的概率来做到这一点：生成 i 作为数字作为 n 个数字序列中的第一个数字的概率，总和为 s 是 F(si, n-1, 0, 6)/F(s, n, 0, 6)。这将保证我们在可能性列表上生成一个 uniform 列表，如果我们每次都绘制一个统一的数字，那么它将不会匹配与给定值匹配的整个值列表上的统一分布条件：

我们可以递归地做到这一点：

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

【讨论】：

这不是以某种奇怪的方式定义的多项分布吗？

@SeverinPappadeux：可能是这样，我会尝试看看。我的想法是（a）实现应该是直截了当的，所以没有“拒绝”，并且（2）好像你会从统一抽样中抽取，然后拒绝。如果我看一下公式，可能就是这种情况，尽管证明它可能需要一些工作:)

如果您发现有趣的事情，请给我留言，谢谢

我更新了我的答案，相信我们提出了不同的解决方案，请看一下

【参考方案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，你的 high 实际上是 high-1

【参考方案3】：

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

顺理成章

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

显然，特定 x₁...x_k 组合的概率会因分母中的阶乘（当然是模排列）而彼此不同，这与@WillemVanOnsem 方法不同，我相信它们都有相同的出现概率。

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