r - R中的指数分布

Question

我想模拟来自 exp(1) 分布的一些数据，但它们必须 > 0.5。所以我使用了一个 while 循环，但它似乎没有像我想要的那样工作。提前感谢您的回复！

x1<-c()

w<-rexp(1) 

while (length(x1) < 100) {

  if (w > 0.5) {

    x1<- w }

  else {

    w<-rexp(1)

  }

}

Answer 1

1）问题中的代码有这些问题：

• 我们在每次迭代中都需要一个新的随机变量，但只有在if条件为 FALSE时才会生成新的随机变量

• x1被反复覆盖而不是扩展

• 虽然while可以使用repeat似乎更好，因为最后的测试比开始的测试更适合

我们可以像这样解决这个问题：

x1 <- c()
repeat {
  w <- rexp(1)
  if (w > 0.5) {
    x1 <- c(x1, w)
    if (length(x1) == 100) break
  }
}

1a)变化如下。请注意，if如果没有else腿，则条件为 FALSE 的 an 评估为 NULL，因此如果在标记为 ## 的行上条件为 FALSE，则没有任何内容连接到x1.

x1 <- c()
repeat {
  w <- rexp(1)
  x1 <- c(x1, if (w > 0.5) w)  ##
  if (length(x1) == 100) break
}

2)或者，这会生成 200 个指数随机变量，仅保留大于 0.5 的变量。如果生成的数量少于 100，则重复。最后，它从生成的最后一批中获取前 100 个。我们选择了足够大的 200，在大多数运行中只需要一次循环迭代。

repeat {
  r <- rexp(200)
  r <- r[r > 0.5]
  if (length(r) >= 100) break
}
r <- head(r, 100)

替代方案 (2) 实际上比 (1) 或 (1a) 更快，因为它的矢量化程度更高。尽管它比其他解决方案丢弃了更多的指数随机变量。

Answer 2

我建议不要使用while（或任何其他接受/拒绝）循环；而是使用以下方法truncdist：

# Sample 1000 observations from a truncated exponential
library(truncdist);
x <- rtrunc(1000, spec = "exp", a = 0.5);

# Plot
library(ggplot2);
ggplot(data.frame(x = x), aes(x)) + geom_histogram(bins = 50) + xlim(0, 10);

使用逆变换采样来实现采样器也相当简单，可以从截断的指数分布中抽取样本，从而避免在循环中拒绝样本。这将是比任何基于接受/拒绝的采样方法更有效的方法，并且在您的情况下效果特别好，因为存在截断指数 cdf 的封闭形式。有关更多详细信息，请参见例如这篇文章。