r - 加快循环中指数随机变量的重复生成

Question

我正在实现一种算法，作为其中的一部分，我需要生成指数随机变量。不幸的是，我无法真正避免循环，因为每个生成的随机变量都依赖于前一个随机变量，所以我认为矢量化是不可能的。我围绕一代做了一些计算，但瓶颈(目前)是一代。在这一点上，我假设 N 会很大(N >= 1,000,000)。

下面是一些示例代码:

N <- 1e7
#Preallocate
x <- rep(0, times=N)

#Set a starting seed
x[1] <- runif(1)

for(i in 2:N) {
    #Do some calculations
    x[i] <- x[i-1] + rexp(1, x[i-1])  #Bottleneck
    #Do some more calculations
}

我怎样才能加快速度？我试过在 Rcpp 中实现，但在这种情况下它似乎并没有做太多。有没有另一种聪明的方法可以绕过每次迭代中的 rexp() 调用？

Answer 1

我们可以利用if X ~ Exp(λ) then kX ~ Exp(λ/k) (source: Wikipedia)加快代码。这样我们就可以预先使用 rate = 1 进行所有随机抽取，然后只需在循环内划分以适本地缩放它们。

draws = rexp(N, rate = 1)
x <- rep(0, times = N)
x[1] <- runif(1)
for(i in 2:N) {
    #Do some calculations
    x[i] <- x[i-1] + draws[i] / x[i-1] 
    #Do some more calculations
}

具有 N = 1e6 值的微基准测试显示这大约快 14 倍:

N <- 1e6
draws = rexp(N, rate = 1)
x <- rep(0, times = N)
x[1] <- runif(1)

microbenchmark::microbenchmark(
    draw_up_front = {
        draws = rexp(N, rate = 1)
        for (i in 2:N)
            x[i] <- x[i - 1] + draws[i] / x[i - 1]
    },
    draw_one_at_time = {
        for (i in 2:N)
            x[i] <- x[i - 1] + rexp(1, x[i - 1])
    },
    times = 10
)
# Unit: milliseconds
#              expr       min        lq      mean    median        uq       max neval cld
#     draw_up_front  153.9547  156.6552  159.9622  160.1901  161.9803  167.2831    10  a 
#  draw_one_at_time 2207.1997 2212.0460 2280.1265 2236.5197 2332.9913 2478.5104    10   b