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

我建议不要使用 while（或任何其他 accept/reject）循环；而是使用 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);

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

Answer 2

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) 下面是一个变体。请注意，如果没有 else 分支，则条件为 FALSE 的 if 的计算结果为 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) 更快，因为它的矢量化程度更高。尽管它比其他解决方案丢弃了更多的指数随机变量。

R 中的指数分布

Exponential distribution in R

r

exponential-distribution