-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です。次のように続けてください。

  • セットnum <- 100(つまり、pp の定義内の「反復」の数)
  • (pai, alpha)長さの次元が追加された配列としてそれぞれを生成しますnum(つまり、各要素またはスライスは、それらの変数の元の反復と同等になります)。
  • taでは毎回 2 つの行列を操作しているため、残りの行列演算はボトルネックです(ta, arrayta, tare, predarray)。
  • ボトルネックは、行列オペランドの各ペアを同じ配列 (または、同じ次元でない場合はリスト内の要素としての各パリ) の追加の次元で生成 (または結合) し、マージンを超えて操作することで修正できます。 「反復」を表します。

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 に答える