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

Question

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

我也需要生成一组正整数随机变量 {Xi}，这些变量加起来就是给定的和 S，其中每个变量可能都有约束，例如 mi

我知道怎么做，问题是在我的例子中，随机变量之间也可能有约束，Xi说

假设我们甚至看一个非常简单的情况：4个随机变量X1、X2、X3、X4需要加起来等于100并符合约束X1

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

Answer 1

我对我的问题有一个答案，在一系列一般限制下，我所做的是：

< li >对约束进行采样，以评估s，即约束区域。 < li >如果s足够大，则生成随机样本，并丢弃那些不符合我之前答案中所述限制的样本。 < li >否则: < ol > < li >枚举整个单纯形。 < li >应用约束以过滤出约束区域之外的所有元组。 < li >列出过滤后得到的元组。 < li >当要求生成时，我从结果列表中选择统一生成。(注意:这值得我努力，因为我经常被要求生成)

注意：我还必须处理 S 是随机生成的参数 （m

Answer 2

好的，对于我的实际问题，我有这个解决方案，我通过将0和排序的随机元组数组连接起来，最后是1，然后按照我在最初问题中提到的答案中的建议，取它们的差，来生成9000个3个随机变量的三元组。

然后我简单地过滤掉那些与我的约束不匹配的并绘制它们。

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/a次）。

如果有人能回答这最后一个问题，我将不胜感激(将选择的答案改为他的)。

Answer 3

对于4个随机变量X1、X2、X3、X4，它们需要加起来为100并符合约束X1

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

代码，Windows 10 x64，Python 3.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

你看到区别了吗？只要应用任何类型的过滤器，就可以更改分布。没有什么是统一的了