生成两个具有固定均值和 sd 的随机数序列，具有排序约束
Generate two sequences of random numbers with fixed mean and sd, with ordering constraint
我想生成具有以下约束的（伪）随机正态分布数字的两个序列 S 和 T：
- S 的均值 m_S、标准差 sd_S 和长度 N
- T 的均值 m_T、标准差 sd_T 和长度 N
- 对于每个i，T(i) <= S(i)，即序列T中位置i的元素小于等于序列S中位置i对应的元素
理想情况下，我想逐步生成两个序列（一次生成一对或数字）。
我知道如何独立生成 S 或 T（例如 Java mean + stdev * ThreadLocalRandom.current().nextGaussian();- 另请参阅 this question），但我不知道如何生成 S 和 T满足第三个约束。
感谢任何帮助！
按照措辞，这在数学上是不可能的。令 D = S - T。您的约束 T <= S 意味着 S - T = D >= 0。由于 S 和 T 呈正态分布，因此 D 也呈正态分布，因为法线的线性组合是正态的。您不能拥有具有有限下限的正态分布。因此，您无法同时满足 S 和 T 的正态性要求并满足您的约束条件。
您可以构造非正态解，方法是生成具有您选择的任何分布的 T，独立生成具有严格正分布（例如伽玛、威布尔、均匀 (0,?)、截断正态分布、... )，并创建 S = T + D。由于独立随机变量的期望值和方差相加，您可以获得所需的均值和 s.d。通过适当地调整 D 的参数化来为 T 和 S。结果甚至看起来很像钟形，但严格来说是不正常的。
由于独立随机变量的方差是可加的，必须是正的，所以S = T + D只有当S的方差大于T的方差时才有效。更一般的解决方案是生成S和T中的那个具有较小的方差。如果是T，加D得S。如果是S，减去D得T。
既然你在评论中说近似值是可以的，这里有一个例子。假设您希望较小的分布具有 μsmaller = 10 和 σsmaller = 3，而较大的分布具有 μ larger = 15 和 σlarger = 5。那么它们之间的差异应该严格正分布，μdelta = 5 σdelta = 4 (σlarger = sqrt(32 + 4 2) = 5).我为 delta 选择了伽玛分布，参数化为具有所需的均值和标准差。这是在 Python:
import random
alpha = 25.0 / 16.0
beta = 16.0 / 5.0
for _ in range(100000):
smaller = random.gauss(10, 3)
delta = random.gammavariate(alpha, beta)
print(smaller, delta, smaller + delta)
我将结果保存到一个文件中，并将它们导入到 JMP 中。这是我的分析快照：
如您所见，较小和较大具有所需的均值和标准差。您还可以确认所有增量都是正的，因此较大的总是大于较小的。最后，较大的直方图上方的正常 q-q 图显示结果虽然是单峰的且大致呈钟形，但并不正常，因为绘制的点不沿直线落下。
另一个答案提出通过生成单个随机均匀分布并将其用作两个 CDF 的反演输入来匹配两个分布：
q = random()
t = inverseCDF(q, mu_T, sd_T)
s = inverseCDF(q, mu_S, sd_S)
这是一种著名的相关归纳策略，称为“Common Random Numbers”，其中 q 是 相同 分位数，用于通过反演生成两个分布。对于对称分布，例如正态分布，这会在 T 和 S 之间产生 1 的相关性。1 的相关性告诉我们（或者）一个是另一个的线性变换。
事实上，有一种更简单的方法可以对法线完成此操作，而无需进行两次反转。通过您希望的任何机制生成 T 或 S 之一——反转、极坐标法、Ziggurat 法等。然后使用标准变换将其转换为标准法线，然后从那里转换为另一个法线。如果我们让 T 是我们直接生成的法线，那么
S = (σS / σT) * (T - μT ) + μS.
我们希望所有可能的分位数都满足 T <= S 以满足原始问题的目标。那么在什么情况下我们可以有S < T呢？由于 S 是 T 的函数，这意味着
(σT/σT) * (T - μT) + μS < T
经过一些代数运算后变成
T * (σS - σT) / σS < μT - μS * (σT/σS ).
这减少到 3 个案例。
σS = σT：在这种情况下，T被淘汰了，T原本想要的结果<= S 只要μT <= μS.
σS > σT：在这种情况下，T > S 当 T < (μT * σS / (σS - σT) ) - (μS * σT / (σS - σT )).
σS < σT: T > S 当 T > (μT * σS / (σS - σT)) - (μ S * σT / (σS - σT)) 由于 (σS - σT) < 0.[=15 在 #2 的结果中引起符号翻转=]
底线 - 相关归纳方案起作用的唯一情况是两个分布的方差相等。不等方差将导致结果 T > S.
下面的图片可能会给出一些直觉。红色曲线是均值为 0 标准差为 1 的标准正态曲线。绿色曲线是均值为 1 标准差为 2 的正态曲线。我们可以看到，由于绿色曲线较宽，因此存在一些分位数，低于该分位数会产生较小的结果比红色。如果 "lower" 分布 T 具有较大的可变性，则会有一些分位数，高于该分位数会产生更大的结果。
我认为你可以计算这样的序列，满足所有条件并递增地生成两个序列（一次一对或数字）。为此，您必须构建累积分布函数 CDFS(x, mS, sdS)和 CDFT(x, mT, sdT) 并确保对于给定的参数ANY x 它们的差异是非负的：CDFS(x) - CDFT(x) >= 0.
然后你抛出 U(0,1)，对于每个均匀分布的随机数，你使用逆 CDF 方法（在一些伪代码中）计算了两个高斯分布：
r = random()
t = inverseCDF(r, mu_T, sd_T)
s = inverseCDF(r, mu_S, sd_S)
它们将成对并保证是高斯分布并满足条件 1,2,3。
逆 CDF 将通过逆误差函数表示，因此根据您使用的库，您可能已经拥有它。如果不是，在 Wiki 中有一些 erf-1(x) 的近似值可以使用：https://en.wikipedia.org/wiki/Error_function
更新
为了更清楚。如果您按照我的建议进行操作并生成，比如说，具有 N 条记录和 2 列（S，T）的数据框，那么显然单独使用 S 列将是高斯（m_S，sd_S）和 T列将是高斯 (m_T, sd_T)。没有问题。他们不会独立。但是你没有说明这样的条件。如果没问题，那就是解决方案。如果他们必须独立，那是不可能的。所以你要拿定主意你到底想要什么
我想生成具有以下约束的（伪）随机正态分布数字的两个序列 S 和 T：
- S 的均值 m_S、标准差 sd_S 和长度 N
- T 的均值 m_T、标准差 sd_T 和长度 N
- 对于每个i，T(i) <= S(i)，即序列T中位置i的元素小于等于序列S中位置i对应的元素
理想情况下，我想逐步生成两个序列（一次生成一对或数字）。
我知道如何独立生成 S 或 T（例如 Java mean + stdev * ThreadLocalRandom.current().nextGaussian();- 另请参阅 this question），但我不知道如何生成 S 和 T满足第三个约束。
感谢任何帮助！
按照措辞，这在数学上是不可能的。令 D = S - T。您的约束 T <= S 意味着 S - T = D >= 0。由于 S 和 T 呈正态分布，因此 D 也呈正态分布，因为法线的线性组合是正态的。您不能拥有具有有限下限的正态分布。因此，您无法同时满足 S 和 T 的正态性要求并满足您的约束条件。
您可以构造非正态解，方法是生成具有您选择的任何分布的 T，独立生成具有严格正分布（例如伽玛、威布尔、均匀 (0,?)、截断正态分布、... )，并创建 S = T + D。由于独立随机变量的期望值和方差相加，您可以获得所需的均值和 s.d。通过适当地调整 D 的参数化来为 T 和 S。结果甚至看起来很像钟形，但严格来说是不正常的。
由于独立随机变量的方差是可加的，必须是正的，所以S = T + D只有当S的方差大于T的方差时才有效。更一般的解决方案是生成S和T中的那个具有较小的方差。如果是T，加D得S。如果是S，减去D得T。
既然你在评论中说近似值是可以的，这里有一个例子。假设您希望较小的分布具有 μsmaller = 10 和 σsmaller = 3，而较大的分布具有 μ larger = 15 和 σlarger = 5。那么它们之间的差异应该严格正分布，μdelta = 5 σdelta = 4 (σlarger = sqrt(32 + 4 2) = 5).我为 delta 选择了伽玛分布，参数化为具有所需的均值和标准差。这是在 Python:
import random
alpha = 25.0 / 16.0
beta = 16.0 / 5.0
for _ in range(100000):
smaller = random.gauss(10, 3)
delta = random.gammavariate(alpha, beta)
print(smaller, delta, smaller + delta)
我将结果保存到一个文件中，并将它们导入到 JMP 中。这是我的分析快照：
如您所见，较小和较大具有所需的均值和标准差。您还可以确认所有增量都是正的，因此较大的总是大于较小的。最后，较大的直方图上方的正常 q-q 图显示结果虽然是单峰的且大致呈钟形，但并不正常，因为绘制的点不沿直线落下。
另一个答案提出通过生成单个随机均匀分布并将其用作两个 CDF 的反演输入来匹配两个分布：
q = random()
t = inverseCDF(q, mu_T, sd_T)
s = inverseCDF(q, mu_S, sd_S)
这是一种著名的相关归纳策略，称为“Common Random Numbers”，其中 q 是 相同 分位数，用于通过反演生成两个分布。对于对称分布，例如正态分布，这会在 T 和 S 之间产生 1 的相关性。1 的相关性告诉我们（或者）一个是另一个的线性变换。
事实上，有一种更简单的方法可以对法线完成此操作，而无需进行两次反转。通过您希望的任何机制生成 T 或 S 之一——反转、极坐标法、Ziggurat 法等。然后使用标准变换将其转换为标准法线，然后从那里转换为另一个法线。如果我们让 T 是我们直接生成的法线，那么
S = (σS / σT) * (T - μT ) + μS.
我们希望所有可能的分位数都满足 T <= S 以满足原始问题的目标。那么在什么情况下我们可以有S < T呢？由于 S 是 T 的函数，这意味着
(σT/σT) * (T - μT) + μS < T
经过一些代数运算后变成
T * (σS - σT) / σS < μT - μS * (σT/σS ).
这减少到 3 个案例。
σS = σT：在这种情况下，T被淘汰了，T原本想要的结果<= S 只要μT <= μS.
σS > σT：在这种情况下，T > S 当 T < (μT * σS / (σS - σT) ) - (μS * σT / (σS - σT )).
σS < σT: T > S 当 T > (μT * σS / (σS - σT)) - (μ S * σT / (σS - σT)) 由于 (σS - σT) < 0.[=15 在 #2 的结果中引起符号翻转=]
底线 - 相关归纳方案起作用的唯一情况是两个分布的方差相等。不等方差将导致结果 T > S.
下面的图片可能会给出一些直觉。红色曲线是均值为 0 标准差为 1 的标准正态曲线。绿色曲线是均值为 1 标准差为 2 的正态曲线。我们可以看到，由于绿色曲线较宽，因此存在一些分位数，低于该分位数会产生较小的结果比红色。如果 "lower" 分布 T 具有较大的可变性，则会有一些分位数，高于该分位数会产生更大的结果。
我认为你可以计算这样的序列，满足所有条件并递增地生成两个序列（一次一对或数字）。为此，您必须构建累积分布函数 CDFS(x, mS, sdS)和 CDFT(x, mT, sdT) 并确保对于给定的参数ANY x 它们的差异是非负的：CDFS(x) - CDFT(x) >= 0.
然后你抛出 U(0,1)，对于每个均匀分布的随机数，你使用逆 CDF 方法（在一些伪代码中）计算了两个高斯分布：
r = random()
t = inverseCDF(r, mu_T, sd_T)
s = inverseCDF(r, mu_S, sd_S)
它们将成对并保证是高斯分布并满足条件 1,2,3。
逆 CDF 将通过逆误差函数表示，因此根据您使用的库，您可能已经拥有它。如果不是，在 Wiki 中有一些 erf-1(x) 的近似值可以使用：https://en.wikipedia.org/wiki/Error_function
更新
为了更清楚。如果您按照我的建议进行操作并生成，比如说，具有 N 条记录和 2 列（S，T）的数据框，那么显然单独使用 S 列将是高斯（m_S，sd_S）和 T列将是高斯 (m_T, sd_T)。没有问题。他们不会独立。但是你没有说明这样的条件。如果没问题，那就是解决方案。如果他们必须独立，那是不可能的。所以你要拿定主意你到底想要什么