在R中快速生成~10 ^ 9步的随机过程

Question

我有以下任务要执行:

生成公式描述的10 ^ 9步骤:

X(0)=0
X(t+1)=X(t)+Y(t)

哪个Y(t)是具有分布的独立随机变量N(0,1).计算指数t的百分比X(t)为负值.

我尝试了以下代码:

  x<-c(0,0)
  z<-0
  loop<-10^9
  for(i in 2:loop) {
    x[1]<-x[2]
    x[2]<-x[1]+rnorm(1, 0, 1)
    if (x[2]<0) {z<-z+1}
  }

但是,它很慢.我怎样才能加快速度？

Answer 1

通常,对于这样的问题,您可以使用Rcpp包将函数一对一地转换为C++.这应该会带来相当大的加速.

一,R版:

random_sum <- function(loop = 1000) {
  x<-c(0,0)
  z<-0
  for(i in 2:loop) {
    x[1]<-x[2]
    x[2]<-x[1]+rnorm(1, 0, 1)
    if (x[2]<0) {z<-z+1}
  }
  z / loop
}
set.seed(123)
random_sum()
# [1] 0.134

现在的C++版本:

library("Rcpp")
cppFunction("
  double random_sum_cpp(unsigned long loop = 1000) {
    double x1 = 0;
    double x2 = 0;
    double z = 0;
    for (unsigned long i = 2; i < loop; i++) {
      x1 = x2;
      x2 = x1 + Rcpp::rnorm(1)[0];
      if (x2 < 0) z = z+1;
    }
    return z/loop;
  }")

set.seed(123)
random_sum_cpp()
# [1] 0.134

为了完整起见,我们还要考虑提出的矢量化版本:

random_sum_vector <- function(loop = 1000) {
  Y = rnorm(loop)
  sum(cumsum(Y)<0)/loop
}
set.seed(123)
random_sum_vector()
# [1] 0.134

我们看到它为相同的随机种子提供了相同的结果,因此它似乎是一个可行的竞争者.

在基准测试中,C++版本和矢量化版本的表现相似,矢量化版本在C++版本上略显优势:

> microbenchmark(random_sum(100000),
                 random_sum_vector(100000),
                 random_sum_cpp(100000))
Unit: milliseconds
                     expr        min         lq       mean     median         uq       max neval
        random_sum(1e+05) 184.205588 199.859266 209.220232 205.137043 211.026740 274.47615   100
 random_sum_vector(1e+05)   6.320690   6.631704   7.273645   6.799093   7.334733  18.48649   100
    random_sum_cpp(1e+05)   8.950091   9.362303  10.663295   9.956996  11.079513  21.30898   100

但是,矢量化版本会通过内存来降低速度,并会耗尽内存以进行长循环.C++版本几乎不使用内存.

对于10 ^ 9步,C++版本在我的机器上运行大约2分钟(110秒).我没试过R版.根据较短的基准,它可能需要大约7个小时.

> microbenchmark(random_sum_cpp(10^9), times = 1)
Unit: seconds
                 expr      min       lq     mean   median       uq      max neval
 random_sum_cpp(10^9) 110.2182 110.2182 110.2182 110.2182 110.2182 110.2182     1

Answer 2

一种解决方案是采用 @G5W 提出的矢量化，但将其分成更小的块以避免任何内存溢出问题。这为您提供了矢量化解决方案的速度，但通过管理块大小，您可以控制进程使用的内存量。

下面将问题分成 1e+07 块，循环 100 次，总共得到 1e+09。

在第一个块的末尾，记录低于 0 的时间百分比以及结束点。然后，结束点将被输入到下一个块，并记录低于 0 的时间百分比以及新的结束点。

最后，对 100 次运行进行平均，使总时间低于零。while 循环中的调用cat是为了监视进度并查看进度，这可以注释掉。

funky <- function(start, length = 1e+07) {
  Y <- rnorm(length)
  Z <- cumsum(Y)
  c(sum(Z<(-start))/length, (tail(Z, 1) + start))
}

starttime <- Sys.time()
resvect <- vector(mode = "numeric", length = 100)
result <- funky(0)
resvect[1] <- result[1]
i <- 2
while (i < 101) {
  cat(result, "\n")
  result <- funky(result[2])
  resvect[i] <- result[1]
  i <- i + 1
}
mean(resvect)
# [1] 0.1880392
endtime <- Sys.time()
elapsed <- endtime - starttime
elapsed
# Time difference of 1.207566 mins