r - 如何矢量化三重嵌套循环？

Question

我已经完成了对类似问题的搜索，但我对我应该做什么有一个模糊的想法:矢量化所有内容或使用 apply()家庭。但我是 R 编程的初学者，上述两种方法都令人困惑。

这是我的源代码:

x<-rlnorm(100,0,1.6)
j=0
k=0
i=0
h=0
lambda<-rep(0,200)
sum1<-rep(0,200)
constjk=0
wj=0
wk=0
for (h in 1:200)
{
   lambda[h]=2+h/12.5
   N=ceiling(lambda[h]*max(x))
   for (j in 0:N)
   {
      wj=(sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
      for (k in 0:N)
      {
         constjk=dbinom(k, j + k, 0.5)
         wk=(sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100
         sum1[h]=sum1[h]+(lambda[h]/2)*constjk*wk*wj
      }
   }
}

让我解释一下。我想收集 200 个 sum1 值(这是第一个循环)，对于每个 sum1 值，它是 (lambda[h]/2)*constjk*wk*wj 的总和，因此其他两个循环。最乏味的是 N 随 h 变化，所以我不知道如何矢量化 j-loop 和 k-loop。但当然我可以用 lambda<-seq() 向量化 h 循环和 N<-ceiling() ，这是我能做的最好的事情。有没有办法进一步简化代码？

Answer 1

您的代码可以完美地使用 3 个嵌套 sapply 进行验证。调用。对于未经训练的眼睛来说可能有点难以阅读，但它的本质是不是一次向 sum1[h] 添加一个值。我们一次性计算出最内层循环产生的所有项并将它们相加。

尽管此矢量化解决方案比您的三元组更快 for循环，改进并不显着。如果你打算多次使用它，我建议你用 C 或 Fortran 实现它(使用常规 for 循环)，这大大提高了速度。但请注意，它具有很高的时间复杂度，并且会随着 lambda 值的增加而严重扩展。 ，最终达到无论实现如何都无法在合理时间内进行计算的地步。

lambda <- 2 + 1:200/12.5
sum1 <- sapply(lambda, function(l){
    N <- ceiling(l*max(x))
    sum(sapply(0:N, function(j){
        wj <- (sum(x <= (j+1)/l) - sum(x <= j/l))/100
        sum(sapply(0:N, function(k){
            constjk <- dbinom(k, j + k, 0.5)
            wk <- (sum(x <= (k+1)/l) - sum(x <= k/l))/100
            l/2*constjk*wk*wj
        }))
    }))
})

顺便说一句，您不需要预定义像 h 这样的变量, j , k , wj和 wk .特别是因为不是在矢量化时，因为在提供给 sapply 的函数中对它们进行赋值将创建具有相同名称的覆盖局部变量(即忽略您预先定义的那些)。