加速嵌套在R中的for循环中的while循环答案

【问题标题】：Speeding up while loop nested in a for loop in R加速嵌套在R中的for循环中的while循环
【发布时间】：2012-05-08 10:29:40
【问题描述】：

在 R 中，我有以下示例模块，它重复了一个 for 循环 n 次：

function(n){
#inputs - n - number of results required
    #reserve n spaces for results
    r_num_successes <- 1:n

    #start looping n times
    for(i in 1:n){

        #set first uniform "random" deviate equal to 0.05 and number of successes to 0
        current_unif <- 0.05
        num_successes <- 0

        #start while loop that updates current_unif - it runs as long as 
        #current_unif is less than 0.95, increments num_successes each loop
        while(current_unif < 0.95){

            #set current_unif to a uniform random deviate between the
            #existing current_unif and 1
            current_unif <- runif(1,current_unif)
            num_successes <- num_successes + 1
        }

        #set the i-th element of the results vector to that final num_successes
        #generated by the while loop
        r_num_successes[i] <- num_successes
    }

            #output the mean of all the successes
    return(mean(r_num_successes))
}

当n 变大时，它开始变得非常缓慢。有什么好的优化方法吗？

【问题讨论】：

请用文字描述这个函数应该做什么。看起来您只是在以 5% 的成功率从统一样本中抽取n。
@Andrie 显然他正在尝试提出一种测量 runif 分布的方法，值得提交给 thedailywtf :-)
@CarlWitthoft Ahah。我就在那里，认为这是最模糊的随机采样器竞赛的新条目......
@dplanet 我看到您已将 cmets 添加到您的代码中。这只是将 cmets 添加到非常无意义的代码中。请用文字描述你想要做什么。
@CarlWitthoft - 我懒洋洋地+1。然而，这只是我正在制作的函数的一个更简单的例子，它有更多的行，因此需要优化。我已经评论了脚本 - 我觉得很难解释，因为它除了证明我的观点之外没有任何用处。

标签： r

【解决方案1】：

使用纯 R 无法显着提高速度。字节编译会给您带来小的改进，但您需要转向已编译的代码才能显着提高速度。

更新：这是一个 Rcpp 解决方案，仅适用于 Dirk :)

> nCode <- '
+   int N = as<int>(n);
+   std::vector<double> rns;
+ 
+   RNGScope scope;  // Initialize Random number generator
+ 
+   for(int i=0; i<N; i++) {
+     double current_unif = 0.05;
+     double num_successes = 0;
+     while(current_unif < 0.95) {
+       current_unif = ::Rf_runif(current_unif, 1.0);
+       num_successes++;
+     }
+     rns.push_back(num_successes);
+   }
+ 
+   double mean = std::accumulate(rns.begin(), rns.end(), 0.0) / rns.size();
+   return wrap(mean);  // Return to R
+ '
>
> library(inline)
> nFunRcpp <- cxxfunction(signature(n="int"), nCode, plugin="Rcpp")
> library(compiler)
> nFunCmp <- cmpfun(nFun)
> system.time(nFun(1e5))
   user  system elapsed 
  3.100   0.000   3.098 
> system.time(nFunCmp(1e5))
   user  system elapsed 
  2.120   0.000   2.114 
> system.time(nFunRcpp(1e5))
   user  system elapsed 
  0.010   0.000   0.016

【讨论】：

天哪，乔希发布了 Rcpp 解决方案？超越真棒。我只有温顺的建议：使用rbenchmark::benchmark() 进行计时。哦，也许预先分配了向量。但是很好地使用 STL 来计算均值。让我印象深刻:)
@DirkEddelbuettel：谢谢。我的“非常好的使用 STL 来计算平均值”是搜索 SO 的结果。 ;) 我没有理由使用rbenchmark，因为在这种情况下没有竞争。我还从你的博客中刷了一些代码来开始。 :)

【解决方案2】：

为了完整起见，这是我对@JoshuaUlrich 的建议：

R> res <- benchmark(nFun(1e5L), nFunCmp(1e5L), nFunRcpp(1e5L), nFun2Rcpp(1e5L),
+                  columns = c("test", "replications", "elapsed", "relative"),
+                  replications=10,
+                  order="relative")
R> print(res)
               test replications elapsed  relative
4 nFun2Rcpp(100000)           10   0.117   1.00000
3  nFunRcpp(100000)           10   0.122   1.04274
2   nFunCmp(100000)           10  13.845 118.33333
1      nFun(100000)           10  23.212 198.39316
R>

nFun2Rcpp 只需添加一行：

rns.reserve(N);

并将分配更改为

rns[i] = num_successes;

而不是使用.push_back()，这会使内存分配更高效。

编辑原来这是不准确的，反映了 随机化 算法。如果我为每个添加一组 set.seed()，则两个 C++ 版本之间的时间是相同的。这里没有可衡量的收益。

【讨论】：

你说得对，rbenchmark 的比率要好得多。 :)