7

私は同様の問題を検索しましたが、何をすべきかについて漠然とした考えがあります。すべてをベクトル化するか、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、したがって他の2つのループです。最も退屈なのは、Nがhとともに変化することです。そのため、jループとkループをベクトル化する方法がわかりません。lambda<-seq()しかしもちろん、とでhループをベクトル化することはできますがN<-ceiling()、それが私にできる最善の方法です。コードをさらに単純化する方法はありますか？

4

2 に答える 2

5

コードは、3 つのネストされた呼び出しで完全に verctorized できsapplyます。慣れていない人には少し読みにくいかもしれませんが、その本質は、一度に 1 つの値を追加するのではなく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の代入は、同じ名前のオーバーレイされたローカル変数を作成するためです (つまり、事前定義したものは無視されます)。

于 2012-11-09T12:39:47.270 に答える
2

シミュレーションを関数でラップして時間を計りましょう。

sim1 <- function(num=20){
  set.seed(42)
  x<-rlnorm(100,0,1.6)
  j=0
  k=0
  i=0
  h=0
  lambda<-rep(0,num)
  sum1<-rep(0,num)
  constjk=0
  wj=0
  wk=0

  for (h in 1:num)
  {
    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)
      {
        set.seed(42)
        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
      }
    }
  }

  sum1
}

system.time(res1 <- sim1())
#   user  system elapsed 
#    5.4     0.0     5.4

それでは、高速化してみましょう。

sim2 <- function(num=20){
  set.seed(42) #to make it reproducible
  x <- rlnorm(100,0,1.6)

  h <- 1:num
  sum1 <- numeric(num)
  lambda <- 2+1:num/12.5
  N <- ceiling(lambda*max(x))

  #functions for wj and wk
  wjfun <- function(x,j,lambda,h){
    (sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
  }
  wkfun <- function(x,k,lambda,h){
    (sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100
  }

  #function to calculate values of sum1
  fun1 <- function(N,h,x,lambda) {
    sum1 <- 0
    set.seed(42) #to make it reproducible
    #calculate constants using outer
    const <- outer(0:N[h],0:N[h],FUN=function(j,k) dbinom(k, j + k, 0.5))
    wk <- numeric(N[h]+1)
    #loop only once to calculate wk
    for (k in 0:N[h]){
      wk[k+1] <- (sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100 
    }

    for (j in 0:N[h])
    {
      wj <- (sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
      for (k in 0:N[h])
      {
        sum1 <- sum1+(lambda[h]/2)*const[j+1,k+1]*wk[k+1]*wj
      }
    }
    sum1
  }

  for (h in 1:num)
  {
    sum1[h] <- fun1(N,h,x,lambda)
  }  
  sum1
}

system.time(res2 <- sim2())
#user  system elapsed 
#1.25    0.00    1.25 

all.equal(res1,res2)
#[1] TRUE

比較のための @Backlin のコード (20 回の反復) のタイミング:

   user  system elapsed 
   3.30    0.00    3.29

これでも遅すぎて、別の言語を使用できない、または使用したくない場合は、並列化の可能性もあります。私が見る限り、外側のループは恥ずかしいほど並列です。並列化のための便利で簡単なパッケージがいくつかあります。

于 2012-11-09T14:13:21.890 に答える