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

Question

我有以下任务要执行:

Generate 10^9 steps of the process described by the formula:

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

where Y(t) are independent random variables with the distribution N(0,1). Calculate in what percentage of indices t the value of X(t) was negative.

我尝试了以下代码:

  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

但是，矢量化版本会牺牲内存和 will blow up your memory for long loops. 的速度。 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