如何生成 3 个在各自特定范围内但总和为 1 的随机变量

Question

我有 3 个变量：X、Y、Z。我想使用 Python 生成 X、Y、Z 的随机数，其中 X+Y+Z=1。

然而，每个变量的范围不同，例如，X 的范围必须在 0.71 到 0.72 之间，Y 的范围必须在 0.23 到 0.24 之间，Z 的范围必须在 0.05 到 0.06 之间。

变量的精度并不重要，因此您可以生成数字： X= 0.71613464 Y=0.23166405 Z=0.05220131

Answer 1

问题

我们想要选择特定范围内的 3 个坐标。然而，一旦我们选择X和Y，Z就固定了。

我们不能拥有完全的自由度，例如，如果我们选择X = 0.72和Y = 0.24，那么就Z不能存在于所需的范围内。

如果我们随机选择X和Y并计算Z，在许多情况下Z将是不正确的：

size = 10_000

X = np.random.uniform(0.71, 0.72, size=size)
Y = np.random.uniform(0.23, 0.24, size=size)
Z = 1-(X+Y)

问题是：如何解决这个问题？

过采样

这通常不是推荐的方法，但我们可以使用试错法，这里我们需要采样两次以上，以希望有足够的有效点（在正方形上采样，保留一半的点）：

target = 100_000
factor = 2.2  # the exact value depends on the target
              # there is always a risk that we miss a few points

size = int(target*factor)
X = np.random.uniform(0.71, 0.72, size=size)
Y = np.random.uniform(0.23, 0.24, size=size)
Z = 1-(X+Y)
idx = np.arange(size)[(Z>=0.05) & (Z<=0.06)][:target]
X = X[idx]
Y = Y[idx]
Z = Z[idx]

assert np.allclose(X+Y+Z, 1)

使用数学

我们正在寻找的相当于在三角形上生成均匀点。

您可以使用这个公式：

size = 10_000
limits = [(0.71, 0.72), (0.23, 0.24), (0.05, 0.06)]

MIN, MAX = np.array(limits).T
corners = np.repeat(MIN[None], 3, axis=0)
np.fill_diagonal(corners, MAX)

# or setting the corners manually:
# corners = np.array([[0.72, 0.23, 0.05],
#                     [0.71, 0.24, 0.05],
#                     [0.71, 0.23, 0.06]])

r1 = np.sqrt(np.random.random(size=size))
r2 = np.random.random(size=size)

X, Y, Z = (corners[:,None] * np.array([(1-r1), r1*(1-r2), r2*r1])[...,None]
          ).sum(axis=0).T

assert np.allclose(X+Y+Z, 1)

einsum最后一步的替代使用：

X, Y, Z = np.einsum('ik,ij->kj', corners, np.array([(1-r1), r1*(1-r2), r2*r1]))

（黑点是三角形的角）

请注意，要使此方法发挥作用，两个坐标的最小边界 + 第三个坐标的最大边界之和必须等于 1 ( Xmax+Ymin+Zmin == Xmin+Ymax+Zmin == Xmin+Ymin+Zmax == 1)。如果不是这种情况，请忽略 MAX 并使用corners = (1-np.identity(3))*MIN ; np.fill_diagonal(corners, 1-corners.sum(axis=1)).

Answer 2

import random

X, Y, Z = 0.0, 0.0, 0.0

epsilon = 0.01

target = 1

while True:
    
    X = random.uniform(0.71, 0.72)
    Y = random.uniform(0.23, 0.24)
    Z = random.uniform(0.05, 0.06)
    
    if target - epsilon <= (X + Y + Z) <= target + epsilon:
        print(X, Y, Z)

更新：

import random

epsilon = 1000

x_d, x_u = 71 * epsilon, 72 * epsilon
y_d, y_u = 23 * epsilon, 24 * epsilon
z_d, z_u = 5 * epsilon, 6 * epsilon


target = 1

while True:
    
    X = random.randint(x_d, x_u)
    Y = random.randint(y_d, y_u)
    Z = random.randint(z_d, z_u)
    
    if (X + Y + Z) == target*epsilon:
        print(X / epsilon, Y / epsilon, Z / epsilon)