-1

我有以下 R 代码：

pp = function(N,J,K){
    for(i in 1:100){
        pai=runif(J)
        alpha=matrix(rbinom(N*K,1,0.5),nrow=N)
        Q=matrix(rbinom(J*K,1,0.5),nrow=J)
        r=matrix(runif(J*K),nrow=J)
        ta=r^Q
        arrayalpha=array(rep((1-alpha), J),c(N,K,J))
        arrayta=array(rep(ta, N),c(J,K,N))
        arraytap=aperm(arrayta, c(3,2,1))
        tare=arraytap^arrayalpha #ta^re
        arrayprod=apply(tare,c(1,3),prod)
        repai=t(matrix(rep(pai,N),nrow=J))
        predarray=t(arrayprod*repai)
   }
   predarray             
}

> system.time(pp(500,20,5))
   user  system elapsed 
  5.381   0.008  12.468

我怎样才能使它更有效率？谢谢您的帮助。

4

1 回答 1

1

考虑到predarray每个循环都被完全替换，您将获得相同的结果而无需迭代（但速度更快）。

不过，说真的，你不想积累 predarray 吗？，比如：

pp = function(N,J,K){
    predarray <- array(numeric(100*J*N), c(100, J, N))
    for(i in 1:100){
        ...
        predarray[i, , ] <- t(arrayprod*repai)}
    predarray             
}

这将返回 100 个 predarray(JxN) 矩阵（每个矩阵作为 (100xJxN) 数组的一个 dim:1 切片）。此外，它应该减少您的时间，因为数组是预先分配的（而不是分配它 100 次，如在您的代码中）。

为了提高效率，您的目标应该是避免显式的 R 循环。最好的方法是显式并行化，因为这个问题看起来“令人尴尬地并行”（即每次迭代都应该完全独立于其他迭代）。

另一种选择是，由于显式 R 循环比隐式 C 循环（用于向量化操作，apply函数族）和包plyr. 你应该继续这样：

  • Set num <- 100（即“迭代”的次数，在pp的定义内）
  • 将每个生(pai, alpha)成为具有附加长度维度的数组num（即每个元素或切片将等效于这些变量的原始迭代）。
  • ta其余的矩阵运算是瓶颈，因为您每次都在两个矩阵上进行运算(ta, arrayta, tare, predarray)。
  • 瓶颈可以通过在同一数组的附加维度中生成（或连接）每对矩阵操作数（或每个 pari 作为列表中的一个元素，当维度不同时），然后在代表“迭代”。

Q 和 r 的示例：

# Q and r in one array of dimensions (J x K x iterations x 2)
# 100 Qs are stored in [,,,1] and 100 rs in [,,,2]
Qr <- array(c(Q=rbinom(J*K*100,1,0.5),
              r=runif(J*K*100)),
            dim=c(J=J, K=K, iter=100, var=2))

# Third dimension is equivalent to each iteration
Qr[,,1,]

# And you operate each Q^r using apply, over each iteration
library(package=plyr)  # plyr is needed for this
ta <- alply(.data=Qr, .margins=3, .fun=function(x) x[,,2]^x[,,1])
于 2013-02-20T05:34:50.040 回答