0

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

模拟器的主要部分是由一个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}
}

这意味着在每次迭代中，我都会更新 and 的值，r并且s在两个结果之一中，填充我的 vector 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其他方式让我重新编写循环以便能够使用所有内核并加快速度？

4

2 回答 2

5

由于您关于转向 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

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

于 2014-11-12T10:35:37.257 回答
1

您正在尝试做的事情，因为每次迭代都依赖于循环的前面步骤，似乎不是可并行化的。您正在更新变量 r 并期望同时运行的其他分支知道它，实际上等待更新发生，
1）不会发生。他们不会等待，他们只会获取 r 的当前值，无论他们运行时的值是什么
2）如果它这​​样做了，它将与没有%dopar%

于 2014-11-11T23:45:39.410 回答