如何从总和为 1 的指数分布生成随机数（概率）

Question

考虑一下我想要的 x 随机数总和为 1 并且分布是指数的。当我使用

x<-c(10,100,1000)

a<-rexp(x[3],rate=1)

a<-a/sum(a)

这会改变分布，对吗？

那么有没有人知道一种方法可以使概率仍然呈指数分布？我知道他们不会再完全独立了。

非常感谢！

Answer 1

是的，归一化改变了分布，事实上，不可能精确地达到你想要的。

---

简单的证明

让某些有限 n 的 X₁, ..., X_n 是要生成其值的随机变量。你有两个要求是

1. X_i~Exp(λ) 对于某些 λ>0 和 i=1,…,n.
2. X₁+…+X_n=1.

虽然这两个单独的要求中的每一个都很容易满足，但不可能同时满足这两个要求。原因是指数分布的probability density function 在 [0,∞) 上是正。这意味着每个 X_i 以正概率获得大于 1 的值，这意味着要求 2 并不总是成立。事实上，它以零概率成立。

---

归一化隐含的概率分布

现在您提出了一种直观的方法，从需求 1 开始并执行归一化 Z_i = X_i / (X₁+ …+X_n) 对于每个 i=1,…,n。然而，很少有分布在诸如加法、乘法，尤其是除法等变换下表现良好，因为随机分母很少易于处理。在这种情况下，我们有额外的复杂性，即 Z_i 的分子和分母是相互依赖的。

尽管如此，Z_i 的确切分布的名称实际上是已知的，它是Dirichlet distribution。要看到这一点，note that X_i~Gamma(1,λ)，其中 λ 充当速率参数。接下来，我们看一下狄利克雷分布的definition：我们从 Y_i~Gamma(α_i, θ) 开始，对于 i=1,…,n 和然后，就像您建议的那样，定义 W_i=Y_i / (Y₁+…+Y_n ）。然后 (W₁,…,W_n)~Dirichlet(α_i,…,α_n)。然而，在要求 1 的情况下，对于每个 i=1,…,n，我们有 α_i=1。因此，您的方法导致 (Z₁,…,Z_n)~Dirichlet(1,…,1)。

然后，您可以使用例如MCMCpack 包来模拟其中的值：

library(MCMCpack)
rdirichlet(1, c(1, 1, 1))
#           [,1]      [,2]       [,3]
# [1,] 0.2088649 0.7444334 0.04670173
sum(rdirichlet(1, c(1, 1, 1)))
# [1] 1

现在查看 Dirichlet(1,...,1) 的 probability density function，您会注意到它实际上是恒定的（当为正时）。因此，在某种程度上，您可能会将其视为多元统一的。如果你想一想它是有道理的（例如，想想 x+y=1, x+y+z=1 上的点）。

然而，多元分布在某种程度上是均匀的，但这并不意味着在边际分布方面有相似之处。事实上，可以show 认为它们是 Beta(1, n-1)。

在 Z_i 被限制为 [0,1]

由于对于某些 λ 值，指数随机变量集中在接近于零的位置，人们可能会错误地认为它们实际上具有有限支持。

X_i~Exp(λ)的累积分布函数为1-exp(-λx)。因此 P(X_i∞ 的极限中为 1，但在这种情况下 X 在分布中收敛到 0。因此，我们不能将非退化指数随机变量限制为 [0,1]。但请注意，对于较大的 λ 固定值，1-exp(-λ) 接近于 1，人们可能会错误地认为 X_i 实际上仅限于 [0,1]。

几个琐碎的演示。首先，Z_i（遵循狄利克雷分布）被限制在 [0,1]。

data <- replicate({
  x <- rexp(5)
  z <- x[1] / sum(x)}, n = 100000)
range(data)
# [1] 1.060492e-06 9.633081e-01
plot(density(data, bw = 0.01))

其次，X~Exp(1) 显然取值大于 1。

x <- rexp(10000)
range(x)
# [1] 7.737341e-05 1.005980e+01
mean(x < 1)
# [1] 0.6391
plot(density(x))

---

按正因子缩放

有多个 cmets 提议使用fact，即指数分布在按正因子缩放时是闭合的，因此如果 X ~ Exp(λ)，则 kX ~ Exp(λ/k)。这当然是对的，但不适用于当前情况。原因是 k = X₁+…+X_n 不是一个常数（意味着对于 X_i 的不同实现，k 是不同的)，因此，kX ~ Exp(λ/k) 不成立。现在，如果我们将 k 视为常数（例如 5），则无法保证 Z_i = X_i / 5 将满足您的要求 2。事实上, 约束的概率为 0。

为了清楚地了解正在发生的事情并且不被@MauritsEvers 的经验“证明”所误导，这里有更多细节。

令 (Ω,F,P) 为概率空间。那么X_i:Ω->R;即，X_i 是一个在 R 中取值 X_i(ω) 的函数，其结果为 ω（将它们想象为 set.seed 值）来自 Ω。现在我们确实有这个性质，对于常数 k，kX_i~Exp(λ/k)。然而，常数意味着无论从 Ω 实现的结果 ω 如何，k 的值总是相同的，就好像 k:Ω->R 是一个常数函数。 @MauritsEvers 建议的是 k = X₁+…+X_n。然而，这被视为一个函数，不是恒定的，并且取决于结果 ω。

一些简单的例子展示了这个逻辑是如何失败的：let k=1/X_i。那么 kX_i=1，这是一个退化的随机变量，而不是一个指数变量。类似地，如果 X~N(0,1)，则 kX=1 而不是 kX~N(0,1/X^2)，这将“遵循”X~N(0,1) 给出 kX 的事实~ N(0,k^2) 表示 常数 k.

---

错误的逻辑

现在，上述错误逻辑的起源可以说是错误处理概率概念 + 直接处理 R 中的模拟值。@MauritsEvers 声称如果我们运行

n <- 3
x <- rexp(n)
k <- sum(x)

那么实现的总和k 可以用作上面提到的常数 k 并期望 kX_i~Exp(?)。如上例所示，对n <- 1 的健全性检查已经表明这种论点存在问题，因为那时x / k 只是1——一个退化的随机变量，而不是一个指数变量。据称k <- sum(x) 是一个有效的选择，因为它是许多已经观察到的实现。这实际上就是这个选择无效的原因。在之前的符号中，我们有 k(ω) = X₁(ω)+…+X_n(ω)，所以 k 不是一个常数函数。

另一种看待它的方式是，如果我们将x 视为某种随机的，那么k 与x 的总和一样随机。现在x 和k 都是数字，实现，但在我们要求 R 打印它们之前，我们都不知道它们的值。常数 k 的定义是我们总是知道它的值，而不管 ω 或 set.seed。

最后，作为一个本科练习，可以考虑查看 kX_i 的 CDF：

P(kX_ii

，因此 kX_i~Exp(λ/k)，正如预期的那样。现在采取n <- 2。在这种情况下，我们正在处理

P(X₁ / (X₁ + X₂)

而且我们再也不能如此轻易地摆脱复杂的分母了。当然，我们可以为 Ω 中的某个固定 ω 定义一个常数 k = X₁(ω)+…+X_n(ω)。但是 Z_i = X_i / (X₁(ω)+…+X_n(ω) ) 不再限于 [0,1] 并且要求 2 再次失败。

---

错误的经验“证明”

最后，有人可能会问，为什么@MauritsEvers 的经验“证明”部分（因为模拟 + 拟合 + 假设检验远非理论证明）声称 Z_i 实际上确实遵循指数分布。

这个“证明”的一个关键要素是取lambda <- 1 和n <- 1000，这是一个相对较大的值。在那种情况下，我们有那个

Z_i = X_i/(X₁+…+X_n) ≈ X_i / n * n / (X₁+…+X_n)。

右手边的第二项，根据大数定律，指向 λ——一个固定数——而第一项紧随我们所知的 Exp(λn)。因此，对于较大的 n，我们得到 Z_i 的 近似 为 λExp(λn)。然而，最初的问题不是关于近似或限制分布。

---

总结

我们可以区分以下三种情况：

1. 小号。 (Z₁, …, Z_n) 遵循 Dirichlet(1,…,1) 分布，边际分布不等于指数分布。用指数近似它们会得到任意糟糕的结果。
2. 大号。 (Z₁, ..., Z_n) 仍然遵循 Dirichlet(1,...,1) 分布，并且边缘分布仍然不等于指数分布。然而，用指数近似它们应该会给出完全有效的实际结果。
3. n->∞ 时的限制情况。随着 n 的增长，每个 Z_i 越来越接近 λExp(λn)。然而，正如我们所看到的，λExp(λn) 趋向于退化的随机变量，它们完全等于 0。

Answer 2

来自?rexp

rexp(n, rate = 1)
   [...]
   n: number of observations. If ‘length(n) > 1’, the length is
      taken to be the number required.

所以

x<-c(10,100,1000)
a<-rexp(x,rate=1)

和

一样

rexp(3, rate = 1)

将其归一化为 1 可确保（指数）概率函数满足（指数）概率密度函数的标准。

---

更新

在与@JuliusVainora 进行了一些晦涩的讨论之后，我将证明a 确实呈指数分布。

1. 让我们重新生成数据：

   x <- c(10, 100, 1000)
   set.seed(2018)
   a <- rexp(x[3], rate=1)
   a <- a / sum(a)
   

   为了重现性，我在这里使用了一个固定的随机种子。

2. 我将拟合贝叶斯指数模型来估计 lambda 基于 a 使用 rstan

   library(rstan)
   stan_code <- "
   data {
       int N;
       real x[N];
   }
   
   parameters {
       real lambda;
   }
   
   model {
       x ~ exponential(lambda);
   }
   "
   
   fit <- stan(
       model_code = stan_code,
       data = list(N = length(a), x = a))
   
   fit
   #Inference for Stan model: b690462e8562075784125cf0e71c81e2.
   #4 chains, each with iter=2000; warmup=1000; thin=1;
   #post-warmup draws per chain=1000, total post-warmup draws=4000.
   #
   #          mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
   #lambda 1000.21    0.80 31.11  941.86  978.74  998.95 1020.84 1062.97  1502    1
   #lp__   5907.27    0.02  0.66 5905.52 5907.09 5907.53 5907.71 5907.75  1907    1
   #
   #Samples were drawn using NUTS(diag_e) at Sun Nov  4 01:09:40 2018.
   #For each parameter, n_eff is a crude measure of effective sample size,
   #and Rhat is the potential scale reduction factor on split chains (at
   #convergence, Rhat=1).
   

3. 我们执行 Kolmogorov-Smirnov 检验来比较 a 的经验分布与指数分布的分布与 lambda 从之前的 Stan 模型估计的分布

   ks.test(a, "pexp", summary(fit)$summary[1, 1])
   #
   #   One-sample Kolmogorov-Smirnov test
   #
   #data:  a
   #D = 0.021828, p-value = 0.7274
   #alternative hypothesis: two-sided
   

   p-值为 0.72，我们未能拒绝从两个不同分布中抽取样本的原假设。

---

更新 2

清除 cmets 中的讨论：

1. straightforward（以及更透明的 IMO）证明指数分布族在按正因子缩放时是封闭的 不必调用整个测度理论机械。

2. 更重要的是，让我们回想一下，任何概率密度函数都定义为

   phi(x) = p(x) * N
   

   在哪里

   N = int p(x) 
   

   积分被p(x) 的样本空间占用，这样

   int phi(x) = 1.
   

   是的，p(x) 在phi 和N 的表达式中都是相同的。重要的部分来了：N 仍然是一个常数，因为我们对整个样本空间求和（积分）。

同样，我们通过（已）抽取的样本的常数总和对从指数分布抽取的样本进行归一化。