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

对于 4 个随机变量 X1，X2，X3，X4，它们需要加起来达到 100 并符合约束 X1 <= 2*X2，可以使用多项式分布

一旦第一个数字的概率足够低，你的 条件几乎总是得到满足，如果没有 - 拒绝并重复。 设计多项式分布的总和等于 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)

KBD公司正如我所说，第一个平均值较低，第二个平均值较高。直方图也略有不同

更新三

好吧，狄利克雷，就这样吧。让我们计算滤波器之前和之后生成器的平均值。法典

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

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

好的，所以我为我的实际问题提供了这个解决方案，我通过将零连接到排序的随机元组数组和最后的 1 来生成 9000 个随机变量的三元组，然后按照我在原始问题中提到的 SO 答案中的建议获取它们的差异。

然后，我只需过滤掉与我的约束不匹配的约束并绘制它们。

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}）。在一般情况下，如果约束将采样区域限制为整个单纯形的非常小的子区域，这可能很有用（因为像这样的重采样意味着要从约束区域 a 采样，您需要从单纯形中采样 1/a 倍）。

如果有人对最后一个问题有答案，我将非常有义务（将所选答案更改为他的答案）。

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

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

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