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

Question

假设我想要 x 个总和为 1 且呈指数分布的随机数。当我使用

x<-c(10,100,1000)

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

a<-a/sum(a)

这会改变分布，对吧？

那么有谁知道概率仍然呈指数分布的方法吗？我知道他们以后不会完全独立了。

非常感谢！

Answer 1

来自 ?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 确实呈指数分布。

让我们重新生成数据：
```
x <- c(10, 100, 1000)
set.seed(2018)
a <- rexp(x[3], rate=1)
a <- a / sum(a)
```
我在这里使用固定的随机种子以实现可重复性。

我将根据 a 使用 rstan

拟合贝叶斯指数模型来估计 lambda

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).

我们执行 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

清理评论中的讨论：

straightforward（IMO 更透明）证明指数分布族在正因子 没有 [=88= 的缩放下是封闭的] 必须调用整个测度论机制。

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

phi(x) = p(x) * N

哪里

N = int p(x)

对 p(x) 的样本 space 进行积分，使得

int phi(x) = 1.

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

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

Answer 2

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

直接证明

令 X₁, …, X_n 对于某些有限 n 是您要生成其值的随机变量。您的两个要求是

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

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

归一化隐含的概率分布

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

然而，Z_i的精确分布的名称实际上是已知的，它是Dirichlet distribution. To see that, note that X_i~Gamma(1,λ), where λ acts as the rate parameter. Next, we look at a 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)=1-exp(-λ) 仅在 λ->∞ 的极限为 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))

按正因子缩放

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

为了清楚地了解正在发生的事情，而不是被@MauritsEvers 的经验 "proofs" 误导，这里有一些更多的细节。

设 (Ω,F,P) 为概率 space。那么X_i:Ω->R;即，X_i 是一个在 R 中取值 X_i(ω) 的函数，其结果为 ω（假设它们为 set.seed 值) 来自 Ω。现在我们确实有这个属性，对于常数 k，kX_i~Exp(λ/k)。然而，常数意味着无论 Ω 的实现结果 ω 为何，k 的值始终相同，就好像 k:Ω->R 是一个常数函数。 @MauritsEvers 提出的是 k = X₁+…+X_n。然而，这被视为一个函数，并不是常数，而是取决于结果 ω。

下面是一些证明此逻辑如何失败的简单示例：令 k=1/X_i。那么kX_i=1，是退化随机变量，不是指数随机变量。类似地，如果 X~N(0,1)，则 kX=1 而不是 kX~N(0,1/X^2)，这将 "follow" 从 X~N(0,1)给出 kX ~ N(0,k^2) 对于 constant 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_i <= x) = P(X_i <= x/k) = 1-exp (-λx/k)

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

P(X₁ / (X₁ + X₂) <= x )

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

错误经验"proofs"

最后，有人可能会问，为什么@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)。但是，最初的问题不是关于近似值或限制分布。

总结

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

小n。 (Z₁, …, Z_n) 服从 Dirichlet(1,…,1) 分布，边际分布不等同于指数分布那些。用指数近似它们会给出任意差的结果。
大号。 (Z₁, …, Z_n) 仍然服从 Dirichlet(1,…,1) 分布，边际分布仍然不等同于指数的。然而，用指数近似值应该给出完全有效的结果以用于实际目的。
当n->∞时的极限情况。随着 n 的增长，每个 Z_i 越来越接近 λExp(λn)。然而，正如我们所见，λExp(λn) 趋向于退化的随机变量恒等于零。

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

how to generate random numbers (probabilities) from exponential distribution that sum up to 1

r

random

exponential-distribution

更新

更新 2