python - 创建受约束的随机数？

Question

清理文本：

如何创建 m=5 的随机数加起来，比如 n=100。但是，第一个随机数是 10 < x1 < 30，第二个随机数是 5 < x2 < 20，第三个随机数是 10 < x3 < 25，等等。所以这五个随机数加起来是 100。如何我可以创建这些受限的五个数字吗？

.

[[

相关问题A1）：创建五个加起来等于100的随机数的标准方法是在[0,100]之间采样四个数字，并将边界0和100相加，然后对这六个数字进行排序[0,x1,x2, x3,x4,100]。我寻求的五个随机数是增量。那是，

100 - x[4] = delta 5
x[4]- x[3] = delta 4
x[3]- x[2] = delta 3
x[2]- x[1] = delta 2
x[1] - 0   = delta 1

这五个增量现在加起来为 100。例如，它们可能是 0、1、2、7、90。下面是一些解决这个问题的代码：

total_sum = 100
n = 5
v = numpy.random.multinomial(total_sum, numpy.ones(n)/n)

]]

.

对于我的问题，我不能允许出现较大的间隔，上面的最大点差是 90-7 = 83，这太宽了。因此，我必须指定更紧密的价差，例如 [10,30]。这意味着最大的随机数是 30，这不允许像 83 这样的大点差。

.

[[

相关问题 A2)：创建五个具有相同边界的数字的部分解决方案，10 < x_i < 30，加起来为 100，如下所示：就像在 A1) 中一样，但将下边界 10 添加到增量中。所以我得到了我想要的五个随机数：

100 - x[4] = delta 5 + 10
x[4]- x[3] = delta 4 + 10
x[3]- x[2] = delta 3 + 10
x[2]- x[1] = delta 2 + 10
x[1] - 0   = delta 1 + 10

基本上，我和 A1) 完全一样，但不是从 0 开始，而是从 10 开始。因此，每个数字都有下界 10，但它们没有上界，它可以很大，太大。如何将上限限制为30？这里的问题是如何限制上边界

]]

.

概括地说，我尝试解决的问题类型如下所示：我需要五个随机数加起来为 100，并且我需要为每个数字分别指定边界，例如第一个随机数为 [10,30]，并且然后 [5,10] 为第二个随机数，[15,35] 为第三个随机数，依此类推。它们加起来必须是 100。

但我正在使用的真实数据有大约 100 个数字 x_i (m=50)，所有这些数字加起来大约为 400,000。对于数字 x_i，范围通常为 [3000,5000]。这些数字并不是很准确，我只是想传达一些关于问题大小的信息。目的是进行 MCMC 模拟，因此需要快速生成这些数字。人们提出了非常优雅的解决方案，它们确实有效，但是它们花费的时间太长，所以我不能使用它们。问题仍然没有解决。理想情况下，我想要一个 O(m) 解决方案和 O(1) 内存解决方案。

这个问题不应该是NP难的，感觉不像。应该有一个多项式时间解，对吧？

Answer 1

假设您需要 [10,30] 中的 n_1、[20,40] 中的 n_2、[30,50] 中的 n_3 和 n1+n2+n3=90

如果您需要每个可能的三元组 (n_1, n_2, n_3) 具有相同的可能性，那将是困难的。(20, n_2, n_3) 形式的三元组的个数大于三元组的个数 (10, n_2, n_3)，所以不能一概而论地选择 n_1。

令人难以置信的缓慢但准确的方法是生成正确范围内的所有 5 个随机数，如果总和不正确，则拒绝整个组。

. . . 啊哈！

我找到了一种有效地参数化选择的方法。不过，首先，为简单起见，请注意下限的总和是可能的最小总和。如果从目标数中减去下限的总和，然后从每个生成的数中减去下限，则会出现每个数在区间 [0, max_k-min_k] 内的问题。这简化了数学和数组（列表）处理。让 n_k 是基于 0 的选择，其中 0<=n_k<=max_k-min_k。

总和的顺序是按字典顺序排列的，所有总和首先从 n_1=0（如果有）开始，然后是 n_1==1 总和，依此类推。总和在每个组中按 n_2 排序，然后按 n_3 排序，依此类推。如果您知道有多少总和添加到目标（称为 T），以及有多少总和以 n_1=0、1、2、... 开头，那么您可以在该列表中找到总和数 S 的起始数 n1。然后您可以将问题简化为添加 n_2+n_3+... 以得到 T-n_1，找到总和数 S -（从小于 n_1 的数开始的数原始总和）。

让 pulse(n) 是 n+1 个的列表： (n+1)*[1] 在 Python 术语中。令 max_k,min_k 为第 k 个选择的限制，m_k = max_k-min_k 为基于 0 的选择的上限。然后从第一个数字的选择有 1+m_1 个不同的“和”，而 pulse(m_k) 给出了分布：1 是使每个和从 0 到 m_1。对于前两个选择，有 m_1+m_+1 个不同的和。事实证明，pulse(m_1) 与 pulse(m_2) 的卷积给出了分布。

是时候停止一些代码了：

    def pulse(width, value=1):
        ''' Returns a vector of (width+1) integer ones. '''
        return (width+1)*[value]

    def stepconv(vector, width):
        ''' Computes the discrete convolution of vector with a "unit"
            pulse of given width.

            Formula: result[i] = Sum[j=0 to width] 1*vector[i-j]
            Where 0 <= i <= len(vector)+width-1, and the "1*" is the value
            of the implied unit pulse function: pulse[j] = 1 for 0<=j<=width.
        '''
        result = width*[0] + vector;
        for i in range(len(vector)):
            result[i] = sum(result[i:i+width+1])
        for i in range(len(vector), len(result)):
            result[i] = sum(result[i:])
        return result

这是专门为仅使用“脉冲”数组进行卷积而编码的，因此卷积中的每个线性组合都只是一个和。

这些仅在最终类解决方案的构造函数中使用：

class ConstrainedRandom(object):
    def __init__(self, ranges=None, target=None, seed=None):
        self._rand = random.Random(seed)
        if ranges != None: self.setrange(ranges)
        if target != None: self.settarget(target)

    def setrange(self, ranges):
        self._ranges = ranges
        self._nranges = len(self._ranges)
        self._nmin, self._nmax = zip(*self._ranges)
        self._minsum = sum(self._nmin)
        self._maxsum = sum(self._nmax)
        self._zmax = [y-x for x,y in self._ranges]
        self._rconv = self._nranges * [None]
        self._rconv[-1] = pulse(self._zmax[-1])
        for k in range(self._nranges-1, 0, -1):
            self._rconv[k-1] = stepconv(self._rconv[k], self._zmax[k-1])

    def settarget(self, target):
        self._target = target

    def next(self, target=None):
        k = target if target != None else self._target
        k = k - self._minsum;
        N = self._rconv[0][k]
        seq = self._rand.randint(0,N-1)
        result = self._nranges*[0]
        for i in range(len(result)-1):
            cv = self._rconv[i+1]
            r_i = 0
            while k >= len(cv):
                r_i += 1
                k -= 1
            while cv[k] <= seq:
                seq -= cv[k]
                r_i += 1
                k -= 1
            result[i] = r_i
        result[-1] = k # t
        return [x+y for x,y in zip(result, self._nmin)]
    
    # end clss ConstrainedRandom

将其用于：

ranges = [(low, high), (low, high), ...]
cr = ConstrainedRandom(ranges, target)
seq = cr.next();
print(seq)
assert sum(seq)==target

seq = cr.next(); # get then get the next one.

...ETC。该类可能会被削减一点，但主要的空间开销在 _rconv 列表中，该列表具有存储的卷积。对于 O(NT) 存储，这大约是 N*T/2。

卷积仅使用范围，在相同约束下生成大量随机数，表构建时间“摊销”为零。就 _rconv 列表中的索引数量而言，.next() 的时间复杂度平均约为 T/2 和 O(T)。

---

要了解算法的工作原理，假设有 3 个从零开始的选择序列，最大值为 (5,7,3)，从 0 开始的目标 T=10。在空闲会话中定义或导入脉冲和 stepconv 函数，然后：

>>> pulse(5)
[1, 1, 1, 1, 1, 1]
>>> K1 = pulse (5)
>>> K2 = stepconv(K1, 7)
>>> K3 = stepconv(K2, 3)
>>> K1
[1, 1, 1, 1, 1, 1]
>>> K2
[1, 2, 3, 4, 5, 6, 6, 6, 5, 4, 3, 2, 1]
>>> K3
[1, 3, 6, 10, 14, 18, 21, 23, 23, 21, 18, 14, 10, 6, 3, 1]
>>> K3[10]
18
>>> sum(K3)
192
>>> (5+1)*(7+1)*(3+1)
192

K3[i] 显示了不同选择 n_1、n_2、n_3 的数量，使得 0 <= n_k <= m_k 和 Σ n_k = i。当应用于其中两个列表时，让 * 表示卷积。然后 pulse(m_2)*pulse(m_3) 给出 n_2 和 n_3 之和的分布：

>>> R23 = stepconv(pulse(7),3)
>>> R23
[1, 2, 3, 4, 4, 4, 4, 4, 3, 2, 1]
>>> len(R23)
11

从 0 到 T=10 的每个值都是（几乎）可能的，因此第一个数字可以有任何选择，并且有 R23[T-n_1] 个可能的三元组添加到以 N1 开头的 T=10。因此，一旦您发现有 18 个可能的和加到 10，请生成一个随机数 S = randint(18) 并通过 R23[T:T-m_1-1:-1] 数组倒数：

>>> R23[10:10-5-1:-1]
[1, 2, 3, 4, 4, 4]
>>> sum(R23[10:10-5-1:-1])
18

请注意，该列表的总和是上面 K3[10] 中计算的总和。健全性检查。无论如何，如果 S==9 是随机选择，那么求在不超过 S 的情况下，可以对该数组的多少个前导项求和。这就是 n_1 的值。在这种情况下 1+2+3 <= S 但 1+2+3+4 > S，所以 n_1 为 3。

如上所述，您可以减少问题以找到 n_2。最终数字（本例中为 n_3）将被唯一确定。

Answer 2

首先，定义你的范围：

ranges = [range(11,30), range(6,20), range(11,25)]

然后列举所有可能性：

def constrainedRandoms(ranges):
    answer = []
    for vector in itertools.product(*ranges):
        if sum(ranges) == 100:
            answer.append(vector)
    return answer

等效单线：

answer = [v for v in itertools.product(*ranges) if sum(v)==100]

然后从中选择一个随机元素：

myChoice = random.choice(answer)

编辑：

笛卡尔积（用 计算itertools.product）本身不会占用太多 RAM（这是因为itertools.product返回了一个使用 O(1) 空间的生成器，但正如您所指出的那样需要很多时间）。只有计算列表 ( answer) 才可以。坏消息是，为了使用random.choice，您确实需要一个列表。如果您真的不想使用列表，那么您可能需要提出一个衰减概率函数。这是一个非常简单的概率函数，您可以使用：

def constrainedRandoms(ranges):
    choices = (v for v in itertools.product(*ranges) if sum(v)==100) # note the parentheses. This is now a generator as well

    prob = 0.5
    decayFactor = 0.97 # set this parameter yourself, to better suit your needs
    try:
        answer = next(choices)
    except StopIteration:
        answer = None
    for choice in choices:
        answer = choice
        if random.uniform(0,1) >= prob:
            return answer
        prob *= decayFactor
    return answer

衰减概率允许您增加选择适合您的约束的下一个向量的概率。这样想：你有一堆约束。很有可能，您将拥有相对较少数量的满足这些约束的向量。这意味着每次您决定不使用某个向量时，存在另一个此类向量的概率就会降低。因此，随着时间的推移，您需要逐渐将当前向量返回为“随机选择的向量”。当然，仍然可以在不返回向量的情况下遍历整个循环结构。
请注意，这种衰减概率的想法允许您不必遍历所有候选向量。随着时间的推移，代码返回正在考虑的当前向量的概率增加，从而降低了它继续并考虑其他向量的概率。因此，这个想法也有助于减轻您对运行时的担忧（虽然，我很尴尬地补充，不是很多）

Answer 3

尝试这个：

import random
a = random.randint(10, 30)
b = random.randint(5, 20)
c = random.randint(10, 25)
d = random.randint(5, 15)
e = 100 - (a+b+c+d)

编辑：

以下是生成n随机数列表、给定n范围约束和所需目标总和的方法：

def generate(constraints, limit):
    ans = [random.choice(c) for c in constraints]
    return ans if sum(ans) == limit else generate(constraints, limit)

这是你将如何使用它：

generate([range(10,31),range(5,21),range(10,26),range(5,16),range(10,26)], 100)
=> [25, 20, 25, 12, 18]

但请注意：如果约束不能保证最终达到总和，则会出现无限循环和堆栈溢出错误，例如：

generate([range(1,11), range(10, 21)], 100)
=> RuntimeError: maximum recursion depth exceeded while calling a Python object

Answer 4

这是一个通用的解决方案：

import random
def constrained_rndms(constraints, total):
    result = []
    for x, y in constraints:
        result.append(random.randint(x,y))
    result.append(total - sum(result))
    return result

s = constrained_rndms([(10,20),(5,20),(10,25),(5,15)],100) # -- [19, 5, 16, 15, 45]
sum(s) # -- 100

Answer 5

可以使用两个跨度、四个跨度、八个跨度等（其中跨度是包括其端点的整数范围）来计算使每个可能的总数的方法的数量。将这些数字列在表格中，您可以向后工作以获取样本。例如，假设有 16 个跨度，每个跨度都包括数字 1 到 9。人们发现有 w = 202752772954792 种方法可以使t = 100. 在 1 到 w 的范围内选择一个随机数 r。搜索（或查找）以找到一个数字 J，使得 J>j 的leftways(t-j)*rightways(j))总和超过 r，而 J>j+1 的总和不超过，其中leftways(i)是使用前 8 个跨度生成 i 的方法的数量，以及rightways(j)是使用最后 8 个跨度生成 j 的方法数。使用 i = tj 与前 8 个跨度和 j 与最后 8 个等进行递归。只要只有一种方法可以得出所需的总数，就会出现基本情况。

下面的代码可以通过按宽度降序排序跨度（即首先列出最宽的跨度）并删除一些交换来修改以更有效地运行。请注意，如果 span 不是按宽度降序排列的，则结果向量的顺序将与 spans 不同。此外，可能有可能用二分搜索替换线性for j ...搜索，findTarget但我还没有弄清楚如何在不使用更多空间的情况下这样做。

通过使用对象来存储树结构，而不仅仅是元组，代码可以变得更清晰、更清晰。

使用的空间可能是O(s²·(lg m))如果有 m 个跨度的最大值总计达到 s。所用时间用于O(s²·(lg m))产品总和的初始制表，并且可能O(m+(lg m)·(s/m))用于O(m+(lg m)·s)每个样本。（我没有正确分析空间和时间要求。）在 2GHz 机器上，如果有 16 个相同的跨度 1...10，所示代码每秒产生大约 8000 个样本；如果有 32 个相同的跨度 1...3，则每秒大约 5000 个样本；如果有 32 个相同的跨度 1...30，则每秒大约 2000 个样本。代码后显示了各种目标总数和跨度集的一些示例输出。

from random import randrange
def randx(hi):   # Return a random integer from 1 to hi
    return 1+randrange(hi) if hi>0 else 0

# Compute and return c with each cell k set equal to 
# sum of products a[k-j] * b[j], taken over all relevant j
def sumprods(lt, rt):
    a, b = lt[0], rt[0]
    (la,ma,aa), (lb,mb,bb) = a, b
    if ma-la < mb-lb:           # Swap so |A| >= |B|
        la, ma, aa, lb, mb, bb = lb, mb, bb, la, ma, aa
    lc, mc = la+lb, ma+mb
    counts = [0]*(mc+1)
    for k in range(lc,mc+1):
        for j in range(max(lb,k-ma), min(mb,k-la)+1):
            counts[k] += aa[k-j] * bb[j]
    return (lc, mc, counts)

def maketree(v):
    lv = len(v)
    if lv<2:
        return (v[0], None, None)
    ltree = maketree(v[:lv/2])
    rtree = maketree(v[lv/2:])
    return (sumprods(ltree, rtree), ltree, rtree)


def findTarget(tototal, target, tree):
    global result
    lt, rt = tree[1], tree[2]
    # Put smaller-range tree second
    if lt[0][1]-lt[0][0] < rt[0][1]-rt[0][0]: lt, rt = rt, lt
    (la,ma,aa), (lb,mb,bb) = lt[0], rt[0]
    lc, mc = la+lb, ma+mb
    count = 0
    for j in range(max(lb,tototal-ma), min(mb,tototal-la)+1):
        i = tototal-j
        count += aa[i] * bb[j]
        if count >= target: break
    # Suppose that any way of getting i in left tree is ok
    if lt[1]: findTarget(i, randx(aa[i]), lt)
    else: result += [i]
    # Suppose that any way of getting j in right tree is ok
    if rt[1]: findTarget(j, randx(bb[j]), rt)
    else: result += [j]

spans, ttotal, tries = [(1,6), (5,11), (2,9), (3,7), (4,9), (4,12), (5,8),
                        (3,5), (2,9), (3,11), (3,9), (4,5), (5,9), (6,13),
                        (7,8), (4,11)],  100, 10

spans, ttotal, tries = [(10*i,10*i+9) for i in range(16)], 1300, 10000

spans, ttotal, tries = [(1,3) for i in range(32)],  64, 10000

spans, ttotal, tries = [(1,10) for i in range(16)], 100, 10

print 'spans=', spans
vals = [(p, q, [int(i>=p) for i in range(q+1)]) for (p,q) in spans]
tree = maketree(vals)
nways = tree[0][2][ttotal]
print 'nways({}) = {}'.format(ttotal, nways)
for i in range(1,tries):
    result = []
    findTarget(ttotal, randx(nways), tree)
    print '{:2}: {}  {}'.format(i, sum(result), result)

在下面显示的输出样本中，带有冒号的行包含样本编号、样本总数和样本值数组。其他行显示了跨度集和获得所需总数的方式的数量。

spans= [(1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10), (1, 10)]
nways(100) = 202752772954792
 1: 100  [8, 9, 1, 2, 8, 1, 10, 6, 6, 3, 6, 10, 6, 10, 10, 4]
 2: 100  [2, 6, 10, 3, 1, 10, 9, 5, 10, 6, 2, 10, 9, 7, 4, 6]
 3: 100  [1, 6, 5, 1, 9, 10, 10, 7, 10, 2, 8, 9, 10, 9, 2, 1]
 4: 100  [10, 7, 6, 3, 7, 8, 6, 5, 7, 7, 5, 1, 9, 6, 9, 4]
 5: 100  [7, 1, 10, 5, 5, 4, 9, 5, 3, 9, 2, 8, 6, 8, 10, 8]
    spans= [(1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3), (1, 3)]
nways(64) = 159114492071763
 1: 64  [2, 2, 1, 3, 1, 2, 2, 2, 1, 2, 3, 3, 3, 2, 2, 1, 2, 3, 1, 3, 1, 3, 2, 1, 2, 3, 2, 2, 1, 2, 2, 2]
 2: 64  [1, 2, 1, 1, 1, 3, 3, 3, 2, 1, 1, 2, 3, 2, 2, 3, 3, 3, 1, 2, 1, 2, 2, 3, 2, 2, 1, 3, 1, 3, 2, 2]
 3: 64  [2, 1, 3, 2, 3, 3, 1, 3, 1, 3, 2, 2, 1, 2, 1, 3, 1, 3, 1, 2, 2, 2, 2, 1, 1, 3, 2, 2, 3, 2, 3, 1]
 4: 64  [2, 3, 3, 2, 3, 2, 1, 3, 2, 2, 1, 2, 1, 1, 3, 2, 2, 3, 3, 1, 1, 2, 2, 1, 1, 2, 3, 3, 2, 1, 1, 3]
 5: 64  [1, 1, 1, 3, 2, 2, 3, 2, 2, 1, 2, 2, 1, 2, 1, 1, 3, 3, 2, 3, 1, 2, 2, 3, 3, 3, 2, 2, 1, 3, 3, 1]
spans= [(0, 9), (10, 19), (20, 29), (30, 39), (40, 49), (50, 59), (60, 69), (70, 79), (80, 89), (90, 99), (100, 109), (110, 119), (120, 129), (130, 139), (140, 149), (150, 159)]
nways(1323) = 5444285920
 1: 1323  [8, 19, 25, 35, 49, 59, 69, 76, 85, 99, 108, 119, 129, 139, 148, 156]
 2: 1323  [8, 16, 29, 39, 48, 59, 69, 77, 88, 95, 109, 119, 129, 138, 147, 153]
 3: 1323  [9, 16, 28, 39, 49, 58, 69, 79, 87, 96, 106, 115, 128, 138, 149, 157]
 4: 1323  [8, 17, 29, 36, 45, 58, 69, 78, 89, 99, 106, 119, 128, 135, 149, 158]
 5: 1323  [9, 16, 27, 34, 48, 57, 69, 79, 88, 99, 109, 119, 128, 139, 144, 158]
    spans= [(1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (
1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (
1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (1, 30), (
1, 30), (1, 30), (1, 30)]
nways(640) = 19144856039395888221416547336829636235610525
 1: 640  [7, 24, 27, 9, 30, 23, 30, 27, 28, 29, 2, 30, 28, 19, 7, 27, 10, 2, 21, 23, 24, 2
7, 24, 16, 29, 8, 13, 23, 2, 19, 27, 25]
 2: 640  [30, 2, 17, 28, 30, 16, 5, 1, 26, 24, 22, 19, 26, 16, 16, 30, 27, 15, 19, 30, 15,
 30, 22, 5, 30, 9, 13, 25, 19, 15, 30, 28]
 3: 640  [2, 24, 1, 23, 20, 5, 30, 22, 24, 19, 22, 9, 28, 29, 5, 24, 14, 30, 24, 16, 26, 2
1, 26, 20, 20, 19, 24, 29, 24, 8, 23, 29]
 4: 640  [7, 20, 16, 24, 22, 14, 28, 28, 26, 8, 21, 9, 22, 24, 28, 19, 5, 13, 9, 24, 25, 2
2, 29, 18, 20, 21, 17, 26, 30, 9, 26, 30]