5

我有以下任务要执行：

生成公式描述的过程的 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}
  }

但是，它非常缓慢。我怎样才能加快速度？

4

5 回答 5

9

一般来说，对于这样的问题，您可以使用 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
于 2017-12-23T16:02:50.970 回答
4

这应该快得多，但是十亿的任何事情都可能需要一段时间。使用较小的长度值（例如 10 ^ 6）对此进行测试可能会很好。

length = 10^9
Y = rnorm(length)
sum(cumsum(Y)<0)/length

编辑

根据@user3666197 的评论，我对此进行了测试，他是正确的。此解决方案适用于较小的数字，但一旦步数变得太大，它就会失败。

我针对 OP 的代码测试了我的“矢量化”版本。当随机游走的长度为 10^8 时，我的代码大约需要 7 秒，而 OP 的代码需要 131 秒（在我的笔记本电脑上）。但是，当我将长度增加到 10^9（根据原始问题）时，我的版本导致大量磁盘交换，我不得不终止该进程。该解决方案在 OP 要求的规模下失败。

于 2017-12-23T15:28:46.597 回答
3

一种解决方案是使用@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
于 2017-12-23T17:34:36.620 回答
0

给定随机源在技术上被构建为其他确定性硬件的能力，以满足对生成流的可重复性的要求以及通过某种伪随机生成器算法“生成”随机性的所有条件，这样的源随机性不容易从纯粹[SERIAL]的-转换为任何形式的“公正”-[CONCURRENT]或真正[PARALLEL]的作案方式。

[SERIAL]也就是说，PRG 步骤是重新定义纯代码执行的任何尝试的中心点（阻塞） 。

这不会改变（非）负值的百分比X(t)，而只是确定对于给定的 PRG 硬件实现，没有更短的方法，而是[SERIAL]产生相互（串行）相关值的纯序列。

展开“慢”循环或准（因为这些值仍然是串行相关的）-矢量化处理（R 语言实现功能利用但几乎是硬件 CPU 指令集级别技巧 - 所以不是语言游戏改变者，而是有点规避一些故意缓慢的代码执行构造函数）是最容易发生的事情。

于 2017-12-23T15:32:57.530 回答
0

使用向量通常会比 for 循环产生更好的性能。非常大的数字（即 10^9）的问题是内存限制。由于您只对负指数的最终百分比感兴趣，因此以下将起作用（在 10^9 步中需要几分钟）。

update_state <- function (curr_state, step_size) {
  n <- min(curr_state$counter, step_size)
  r <- rnorm(min(curr_state$counter, step_size))
  total <- curr_state$cum_sum + cumsum(r)

  list('counter' = curr_state$counter - n,
       'neg_count' = curr_state$neg_count + length(which(total < 0)),
       'cum_sum' = curr_state$cum_sum + sum(r))
}


n <- 10^9
curr_state <- list('counter' = n, 'neg_count' = 0, 'cum_sum' = 0)

step_size <- 10^8
while (curr_state$counter > 0) {
  curr_state <- update_state(curr_state = curr_state, step_size = step_size)
}

print(curr_state)
print(curr_state$neg_count/ n)
于 2017-12-23T17:47:12.740 回答