在 R 中定义指数分布以估计概率答案

【问题标题】：Defining exponential distribution in R to estimate probabilities在 R 中定义指数分布以估计概率
【发布时间】：2017-11-20 22:45:27
【问题描述】：

我有一堆随机变量(X1,....,Xn)，它们是 i.i.d。 Exp(1/2) 表示某个事件的持续时间。所以这个分布的期望值显然是 2，但是我在 R 中定义它时遇到问题。我做了一些研究，发现了一些关于所谓的蒙特卡洛刺激的东西，但我似乎没有找到我正在寻找的东西在里面。

我想估计的一个例子是：假设我们有 10 个随机变量 (X1,..,X10) 分布如上，我们想确定例如概率 P([X1+...+X10<=25])。

谢谢。

【问题讨论】：

【解决方案1】：

你知道 R 中的rexp() 函数吗？通过在 R 控制台中输入 ?rexp 查看文档页面。

对期望概率的蒙特卡洛估计的快速回答：

mean(rowSums(matrix(rexp(1000 * 10, rate = 0.5), 1000, 10)) <= 25)

我生成了 1000 组 10 个指数样本，将它们放入 1000 * 10 矩阵中。我们取行总和并得到一个包含 1000 个条目的向量。 0 到 25 之间的值的比例是对所需概率的经验估计。

谢谢，这很有帮助！我可以将replicate 与此代码一起使用，使其看起来像这样：F <- function(n, B=1000) mean(replicate(B,(rexp(10, rate = 0.5)))) 但我无法输出正确的结果。

replicate 这里也生成了一个矩阵，但它是一个 10 * 1000 的矩阵（而不是我的答案中的 1000* 10 矩阵），所以你现在需要使用colSums。还有，你把n放在哪里了？

正确的函数是

F <- function(n, B=1000) mean(colSums(replicate(B, rexp(10, rate = 0.5))) <= n)

对于您给定示例的非蒙特卡洛方法，请参阅其他答案。指数分布是伽马分布的一个特例，伽马分布具有可加性。

我给你蒙特卡洛方法是因为你在你的问题中命名它，它适用于你的例子。

【讨论】：

【解决方案2】：

在这种情况下，您实际上不需要蒙特卡罗模拟，因为：

如果 Xi ~ Exp(λ) 则和 (X1 + ... + Xk) ~ Erlang(k, λ) 这只是 Gamma(k, 1/λ) （在 (k, θ) 参数化中）或 Gamma(k, λ)（在 (α,β) 参数化中），具有整数形状参数 k。

所以，P([X1+...+X10

pgamma(25, shape=10, rate=0.5)

【讨论】：