如何对三重嵌套循环进行矢量化？

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循环和k循环进行矢量化.但是当然我可以使用lambda<-seq()和矢量化h循环N<-ceiling(),这是我能做的最好的.有没有办法进一步简化代码？

Answer 1

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

虽然这种矢量化解决方案比你的tripple 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将创建具有相同名称的重叠局部变量(即忽略您预先定义的那些变量).