在 R 中定义指数分布以估计概率
Defining exponential distribution in R to estimate probabilities
我有一堆随机变量 (X1,....,Xn)，它们是 i.i.d。 Exp(1/2) 表示某个事件的持续时间。所以这个分布显然有 2 的预期值，但我在 R 中定义它时遇到问题。我做了一些研究并发现了一些关于所谓的蒙特卡洛刺激的东西，但我似乎没有找到我正在寻找的东西因为在里面。
我想估计的一个例子是：假设我们有 10 个随机变量 (X1,..,X10) 分布如上，我们想确定概率 P([X1+...+X10<=25]).
谢谢。
你知道 R 中的 rexp() 函数吗？通过在 R 控制台中键入 ?rexp 查看文档页面。
快速回答您 Monte Carlo 所需概率的估计：
mean(rowSums(matrix(rexp(1000 * 10, rate = 0.5), 1000, 10)) <= 25)
我生成了 1000 组 10 个指数样本，将它们放入 1000 * 10 矩阵中。我们取行总和并得到一个包含 1000 个条目的向量。 0 到 25 之间的值的比例是所需概率的经验估计。
Thanks, this was helpful! Can I use replicate with this code, to make it look like this: F <- function(n, B=1000) mean(replicate(B,(rexp(10, rate = 0.5)))) but I am unable to output the right result.
replicate 这里也生成一个矩阵，但它是一个 10 * 1000 矩阵（与我的回答中的 1000* 10 矩阵相反），所以你现在需要 colSums .还有，你把n放在哪里了？
正确的函数应该是
F <- function(n, B=1000) mean(colSums(replicate(B, rexp(10, rate = 0.5))) <= n)
对于给定示例的非Monte Carlo 方法，请参阅其他答案。指数分布是伽玛分布的特例，后者具有可加性属性.
我给你 Monte Carlo 方法是因为你在你的问题中给它命名，而且它适用于你的例子之外。
在这种情况下您实际上不需要 monte carlo 模拟，因为：
If Xi ~ Exp(λ) then the sum (X1 + ... + Xk) ~ Erlang(k, λ) which is just a Gamma(k, 1/λ) (in (k, θ) parametrization) or Gamma(k, λ) (in (α,β) parametrization) with an integer shape parameter k.
来自维基百科 (https://en.wikipedia.org/wiki/Exponential_distribution#Related_distributions)
因此，P([X1+...+X10<=25]) 可以通过
计算
pgamma(25, shape=10, rate=0.5)
我有一堆随机变量 (X1,....,Xn)，它们是 i.i.d。 Exp(1/2) 表示某个事件的持续时间。所以这个分布显然有 2 的预期值，但我在 R 中定义它时遇到问题。我做了一些研究并发现了一些关于所谓的蒙特卡洛刺激的东西，但我似乎没有找到我正在寻找的东西因为在里面。
我想估计的一个例子是：假设我们有 10 个随机变量 (X1,..,X10) 分布如上，我们想确定概率 P([X1+...+X10<=25]).
谢谢。
你知道 R 中的 rexp() 函数吗？通过在 R 控制台中键入 ?rexp 查看文档页面。
快速回答您 Monte Carlo 所需概率的估计：
mean(rowSums(matrix(rexp(1000 * 10, rate = 0.5), 1000, 10)) <= 25)
我生成了 1000 组 10 个指数样本，将它们放入 1000 * 10 矩阵中。我们取行总和并得到一个包含 1000 个条目的向量。 0 到 25 之间的值的比例是所需概率的经验估计。
Thanks, this was helpful! Can I use
replicatewith this code, to make it look like this:F <- function(n, B=1000) mean(replicate(B,(rexp(10, rate = 0.5))))but I am unable to output the right result.
replicate 这里也生成一个矩阵，但它是一个 10 * 1000 矩阵（与我的回答中的 1000* 10 矩阵相反），所以你现在需要 colSums .还有，你把n放在哪里了？
正确的函数应该是
F <- function(n, B=1000) mean(colSums(replicate(B, rexp(10, rate = 0.5))) <= n)
对于给定示例的非Monte Carlo 方法，请参阅其他答案。指数分布是伽玛分布的特例，后者具有可加性属性.
我给你 Monte Carlo 方法是因为你在你的问题中给它命名，而且它适用于你的例子之外。
在这种情况下您实际上不需要 monte carlo 模拟，因为：
If Xi ~ Exp(λ) then the sum (X1 + ... + Xk) ~ Erlang(k, λ) which is just a Gamma(k, 1/λ) (in (k, θ) parametrization) or Gamma(k, λ) (in (α,β) parametrization) with an integer shape parameter k.
来自维基百科 (https://en.wikipedia.org/wiki/Exponential_distribution#Related_distributions)
因此，P([X1+...+X10<=25]) 可以通过
计算pgamma(25, shape=10, rate=0.5)