r - 在 R 中定义指数分布以估计概率

Question

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

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

谢谢。

score 2 · Accepted Answer

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

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

因此，P([X1+...+X10<=25]) 可以计算为

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

score 1 · Accepted Answer

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

快速回答您对所需概率的蒙特卡洛估计：

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 回答 2