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

Question

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

x<-c(10,100,1000)

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

a<-a/sum(a)

这会改变分布，对吧？

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

非常感谢！

score 2 · Accepted Answer

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

直截了当的证明

令 X ₁ , ..., X _n为一些有限的 n 是您想要生成其值的随机变量。你有两个要求是

X _i ~Exp(λ) 对于某些 λ>0 和 i=1,…,n。
X ₁ +…+X _n =1。

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

归一化隐含的概率分布

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

然而，Z _i的精确分布的名称实际上是已知的，它是狄利克雷分布。要看到这一点，请注意X _i ~Gamma(1,λ)，其中 λ 充当速率参数。接下来，我们看一下狄利克雷分布的定义：我们从 Y _i ~Gamma(α _i , θ) for 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) 的概率密度函数，您会注意到它实际上是常数（当为正时）。因此，在某种程度上，您可能会将其视为多元统一的。如果你想一想它是有道理的（例如，想想 x+y=1，x+y+z=1 上的点）。

然而，多元分布在某种程度上是均匀的，但这并不意味着在边际分布方面有相似之处。事实上，可以证明它们是 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))

按正因子缩放

有多个评论建议使用指数分布在按正因子缩放时闭合的事实，因此如果 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 _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 再次失败。

错误的经验“证明”

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

概括

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

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

score 0 · Accepted Answer

从?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)
```
为了重现性，我在这里使用了固定的随机种子。

我将拟合一个贝叶斯指数模型来估计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).

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

要清除评论中的讨论：

证明指数分布族在按正因子缩放时是封闭的，而无需调用整个测量理论机制，这很简单（而且 IMO 更加透明）。
更重要的是，让我们回想一下，任何概率密度函数都定义为
```
phi(x) = p(x) * N
```
在哪里
```
N = int p(x) 
```
积分被接管的样本空间p(x)使得
```
int phi(x) = 1.
```
是的，在 for和 forp(x)的表达式中都是一样的。这是重要的部分：当我们对整个样本空间求和（积分）时，它仍然是一个常数。phiNN

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

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

2 回答 2

更新

更新 2

Related

Reference