生成总和为给定数字并符合一组一般约束的随机自然数

Question

我有一个应用程序需要类似于 here 中描述的问题。

我也需要生成一组加起来等于给定总和 S 的正整数随机变量 {Xi}，其中每个变量可能具有约束条件，例如 mi<=Xi<=Mi。

我知道该怎么做，问题是在我的情况下，我也可能 随机变量本身之间存在约束 ，比如 Xi<=Fi(Xj) 对于一些给定 Fi（也可以说 Fi 的逆已知），Now，应该如何“正确”生成随机变量？我在这里正确地加了引号，因为我不太确定它在这里意味着什么，只是我希望生成的数字能够涵盖所有可能的情况，并且对每种可能的情况都具有尽可能统一的概率。

假设我们甚至看一个非常简单的案例： 4个随机变量X1,X2,X3,X4需要加起来为100并遵守约束条件X1 <= 2*X2,what将是生成它们的“正确”方式吗？

P.S。我知道这似乎更适合数学溢出，但我也没有找到解决方案。

Answer 1

对于4个随机变量X1,X2,X3,X4需要加起来为100且满足约束条件X1 <= 2*X2，可以使用multinomial distribution

只要第一个数字的概率足够低，你的几乎总是会满足条件，如果不满足 - 拒绝并重复。设计的多项式分布总和等于 100。

代码，Windows10 x64，Python3.8

import numpy as np

def x1x2x3x4(rng):
    while True:
        v = rng.multinomial(100, [0.1, 1/2-0.1, 1/4, 1/4])
        if v[0] <= 2*v[1]:
            return v

    return None

rng = np.random.default_rng()

print(x1x2x3x4(rng))
print(x1x2x3x4(rng))
print(x1x2x3x4(rng))

更新

选择概率的自由度很高。例如，您可以使其他 (##2, 3, 4) 对称。代码

def x1x2x3x4(rng, pfirst = 0.1):
    pother = (1.0 - pfirst)/3.0
    while True:
        v = rng.multinomial(100, [pfirst, pother, pother, pother])
        if v[0] <= 2*v[1]:
            return v

    return None

更新二

如果您开始拒绝组合，那么您会人为地增加一个事件子集的概率并降低另一组事件的概率 - 总和始终为 1。没有办法在您想要的条件下获得统一的概率遇到。下面的代码以等概率运行多项式并计算直方图和平均值。平均值应该正好是 25 (=100/4)，但是一旦您拒绝了一些样本，您就会降低第一个值的平均值并增加第二个值的平均值。差异很小，但不可避免。如果你觉得没问题，那就这样吧。代码

import numpy as np
import matplotlib.pyplot as plt

def x1x2x3x4(rng, summa, pfirst = 0.1):
    pother = (1.0 - pfirst)/3.0
    while True:
        v = rng.multinomial(summa, [pfirst, pother, pother, pother])
        if v[0] <= 2*v[1]:
            return v
    return None

rng = np.random.default_rng()

s = 100
N = 5000000

# histograms
first = np.zeros(s+1)
secnd = np.zeros(s+1)
third = np.zeros(s+1)
forth = np.zeros(s+1)

mfirst = np.float64(0.0)
msecnd = np.float64(0.0)
mthird = np.float64(0.0)
mforth = np.float64(0.0)

for _ in range(0, N): # sampling with equal probabilities
    v = x1x2x3x4(rng, s, 0.25)

    q = v[0]
    mfirst   += np.float64(q)
    first[q] += 1.0

    q = v[1]
    msecnd   += np.float64(q)
    secnd[q] += 1.0

    q = v[2]
    mthird   += np.float64(q)
    third[q] += 1.0

    q = v[3]
    mforth   += np.float64(q)
    forth[q] += 1.0

x = np.arange(0, s+1, dtype=np.int32)

fig, axs = plt.subplots(4)
axs[0].stem(x, first, markerfmt=' ')
axs[1].stem(x, secnd, markerfmt=' ')
axs[2].stem(x, third, markerfmt=' ')
axs[3].stem(x, forth, markerfmt=' ')
plt.show()

print((mfirst/N, msecnd/N, mthird/N, mforth/N))

打印

(24.9267492, 25.0858356, 24.9928602, 24.994555)

注意！正如我所说，第一个平均值较低，第二个平均值较高。直方图也有点不同

更新三

好吧，狄利克雷，就这样吧。让我们计算过滤器前后生成器的平均值。代码

import numpy as np

def generate(n=10000):
    uv = np.hstack([np.zeros([n, 1]),
                    np.sort(np.random.rand(n, 2), axis=1),
                    np.ones([n,1])])
    return np.diff(uv, axis=1)

a = generate(1000000)

print("Original Dirichlet sample means")
print(a.shape)
print(np.mean((a[:, 0] * 100).astype(int)))
print(np.mean((a[:, 1] * 100).astype(int)))
print(np.mean((a[:, 2] * 100).astype(int)))

print("\nFiltered Dirichlet sample means")
q = (a[(a[:,0]<=2*a[:,1]) & (a[:,2]>0.35),:] * 100).astype(int)
print(q.shape)

print(np.mean(q[:, 0]))
print(np.mean(q[:, 1]))
print(np.mean(q[:, 2]))

我有

Original Dirichlet sample means
(1000000, 3)
32.833758
32.791228
32.88054

Filtered Dirichlet sample means
(281428, 3)
13.912784086871243
28.36360987535
56.23109285501087

你看出区别了吗？一旦应用任何类型的过滤器，就会改变分布。没有什么是统一的了

Answer 2

好的，所以我有这个解决方案来解决我的实际问题，我生成 9000 个 3 个随机变量的三元组，方法是将零连接到已排序的随机元组数组，最后连接到一个，然后按照 the answer on SO I mentioned in my original question 中的建议获取它们的差异。

然后我简单地过滤掉那些不符合我的约束条件的并绘制它们。

S = 100

def generate(n=9000):
    uv = np.hstack([np.zeros([n, 1]),
                    np.sort(np.random.rand(n, 2), axis=1),
                    np.ones([n,1])])
    return np.diff(uv, axis=1)

a = generate()

def plotter(a):
    fig = plt.figure(figsize=(10, 10), dpi=100)
    ax = fig.add_subplot(projection='3d')

    surf = ax.scatter(*zip(*a), marker='o', color=a / 100)
    ax.view_init(elev=25., azim=75)
    
    ax.set_xlabel('$A_1$', fontsize='large', fontweight='bold')
    ax.set_ylabel('$A_2$', fontsize='large', fontweight='bold')
    ax.set_zlabel('$A_3$', fontsize='large', fontweight='bold')
    lim = (0, S);
    ax.set_xlim3d(*lim);
    ax.set_ylim3d(*lim);
    ax.set_zlim3d(*lim)
    plt.show()

b = a[(a[:, 0] <= 3.5 * a[:, 1] + 2 * a[:, 2]) &\
      (a[:, 1] >= (a[:, 2])),:] * S
plotter(b.astype(int))

如您所见，分布均匀分布在单纯形的这些任意限制上，但我仍然不确定我是否可以放弃不遵守约束的样本（以某种方式将约束处理成生成过程？我现在几乎可以肯定它不能为一般的 {Fi} 完成）。这在一般情况下可能很有用，如果您的约束将采样区域限制为整个单纯形的一个非常小的子区域（因为像这样重新采样意味着要从受限区域采样，您需要从单纯形中采样 1/一次）。

如果有人对最后一个问题有答案，我将不胜感激（会将所选答案更改为他的答案）。

Answer 3

我有一个问题的答案，在一般的约束条件下，我所做的是：

对约束进行采样以评估约束区域 s。
如果 s 足够大，则生成随机样本并丢弃那些不符合我之前回答中描述的约束的样本。
否则：
1. 枚举整个单纯形。
2. 应用约束以过滤掉约束区域外的所有元组。
3. 列出过滤后的元组。
4. 要求生成，我是从这个结果列表中统一选择生成的。（注意：这值得我努力只是因为我经常被要求生成）
这两种策略的组合应涵盖大多数情况。

注意：我还必须处理 S 是随机生成的参数 (m < S < M) 的情况，在这种情况下，我只是将其视为另一个约束在 m 和 M 之间的随机变量，并将它与其余变量并按照我之前的描述进行处理。

生成总和为给定数字并符合一组一般约束的随机自然数

Generate random natural numbers that sum to a given number and comply to a set of general constraints

language-agnostic

random

math

constraints