%dopar%或加速顺序随机计算的替代方法

Question

我写了一个随机过程模拟器,但我想加快它,因为它很慢.

模拟器的主要部分是一个for循环,我想重写为foreach%%dopar%.

我试过用简化的循环这样做,但我遇到了一些问题.假设我的for循环看起来像这样

library(foreach)

r=0
t<-rep(0,500)
for(n in 1:500){
    s<-1/2+r
    u<-runif(1, min = 0, max = 1)
    if(u<s){
        t[n]<-u
        r<-r+0.001
    }else{r<-r-0.001}
}

这意味着在每次迭代时我都会更新和的值,r并且s在两个结果中的一个中填充我的向量t.我已经尝试了几种不同的方式将它重新编写为foreach循环,但似乎每次迭代我的值都没有得到更新,我得到一些非常奇怪的结果.我尝试过使用return但似乎没有用!

这是我想出的一个例子.

rr=0
tt<-foreach(i=1:500, .combine=c) %dopar% {
    ss<-1/2+rr
    uu<-runif(1, min = 0, max = 1)
    if(uu<=ss){
        return(uu)
        rr<-rr+0.001
    }else{
        return(0)
        rr<-rr-0.001}
}

如果不可能使用foreach其他方式让我重新编写循环,以便能够使用所有内核并加快速度？

Answer 1

由于你的评论,关于转向C,是令人鼓舞的,并且 - 最近 - 证明这不是一项艰巨的任务(特别是对于此类操作)并且值得研究,这里是两个接受一些示例函数的示例函数的比较迭代并执行循环的步骤:

ffR = function(n) 
{
   r = 0
   t = rep(0, n)
   for(i in 1:n) {
       s = 1/2 + r
       u = runif(1)
       if(u < s) {
           t[i] = u
           r = r + 0.001
       } else r = r - 0.001
   }

   return(t)
}


ffC = inline::cfunction(sig = c(R_n = "integer"), body = '
    int n = INTEGER(AS_INTEGER(R_n))[0];

    SEXP ans;
    PROTECT(ans = allocVector(REALSXP, n));

    double r = 0.0, s, u, *pans = REAL(ans);

    GetRNGstate();

    for(int i = 0; i < n; i++) {
        s = 0.5 + r;
        u = runif(0.0, 1.0);

        if(u < s) {
            pans[i] = u;
            r += 0.001;
        } else {
            pans[i] = 0.0;
            r -= 0.001;
        }
    }

    PutRNGstate();

    UNPROTECT(1);
    return(ans);
', includes = "#include <Rmath.h>")

比较结果:

set.seed(007); ffR(5)
#[1] 0.00000000 0.39774545 0.11569778 0.06974868 0.24374939
set.seed(007); ffC(5)
#[1] 0.00000000 0.39774545 0.11569778 0.06974868 0.24374939

速度比较:

microbenchmark::microbenchmark(ffR(1e5), ffC(1e5), times = 20)
#Unit: milliseconds
#       expr        min         lq     median         uq        max neval
# ffR(1e+05) 497.524808 519.692781 537.427332 668.875402 692.598785    20
# ffC(1e+05)   2.916289   3.019473   3.133967   3.445257   4.076541    20

为了完整起见:

set.seed(101); ans1 = ffR(1e5)
set.seed(101); ans2 = ffC(1e5)
all.equal(ans1, ans2)
#[1] TRUE

希望这些都可以在某种程度上有所帮助.