2

R で for ループの中に for ループを入れ子にした次のコードを作成しました。 Power を計算するシミュレーションです。R は for ループを実行するのに適していないことを読みましたが、これを少し速く実行するために適用できる効率があるかどうか疑問に思っていました。私はRだけでなく、あらゆる種類のプログラミングにもかなり慣れていません。現在、私が見ている実行時間は次のとおりです。

m=10 0.17 秒

m=100 3.95 秒

m=1000 246.26 秒

m=2000 私は 1003.55 秒を取得します

サンプリングする回数 m を 100K 以上に設定したいと思っていましたが、これを 10K に設定することさえ恐れています

コードは次のとおりです。

m = 1000                        # number of times we are going to  take samples
popmean=120                     # set population mean at 120
popvar=225                      # set known/established population 
variance at 225
newvar=144                      # variance of new methodology 
alpha=.01                       # set alpha
teststatvect = matrix(nrow=m,ncol=1)    # empty vector to populate with test statistics
power = matrix(nrow=200,ncol=1)     # empty vector to populate with power

system.time(                    # not needed - using to gauge how long this takes
    for (n in 1:length(power))          # begin for loop for different sample sizes
      for(i in 1:m){                # begin for loop to take "m" samples
      y=rnorm(n,popmean,sqrt(newvar))   # sample of size n with mean 120 and var=144
      ts=sum((y-popmean)^2/popvar)      # calculate test statistic for each sample
      teststatvect[i]=ts            # loop and populate the vector to hold test statistics
      vecpvals=pchisq(teststatvect,n)   # calculate the pval of each statistic
      power[n]=length(which(vecpvals<=alpha))/length(vecpvals) # loop to populate      power vector. Power is the proportion lessthan ot equal to alpha
        }
   }
 )
4

2 に答える 2

3

コードを少し再編成し、内側のループを取り除きました。

  • 乱数の1つの長いベクトルをサンプリングする（そしてそれを行列に折りたたむ）ことは、短いベクトルを繰り返しサンプリングするよりもはるかに高速です（replicate別の回答で示唆されているように、読みやすさには優れていますが、この場合、乱数をブロック）
  • colSumsforループ内で合計したり、 を使用したりするよりも高速ですapply。
  • それは単なる砂糖です（つまり、実際にはそれ以上効率的ではありません）がmean(pvals<=alpha)、代わりに使用できますsum(pvals<=alpha)/length(alpha)
  • 指定された一連のパラメーター (サンプル サイズを含む) の検出力を返す関数を定義し、sapplyサイズのベクトルを範囲指定するために使用しました (forループよりも高速ではありませんが、よりクリーンで、一般化するのが簡単かもしれません)。

コード：

powfun <- function(ssize=100,
                   m=1000,      ## samples per trial
                   popmean=120, ## pop mean
                   popvar=225,  ## known/established pop variance
                   newvar=144,  ## variance of new methodology
                   alpha=0.01,
                   sampchisq=FALSE)  ## sample directly from chi-squared distrib?
{
    if (!sampchisq) {
      ymat <- matrix(rnorm(ssize*m,popmean,sd=sqrt(newvar)),ncol=m)
      ts <- colSums((ymat-popmean)^2/popvar)          ## test statistic
    } else {
      ts <- rchisq(m,df=ssize)*newvar/popvar
    }
    pvals <- pchisq(ts,df=ssize)                    ## pval
    mean(pvals<=alpha)                              ## power
}

サンプルサイズのすべての整数値のパワーが本当に必要ですか、それともより広い間隔のサンプルでも問題ありませんか(正確な値が必要な場合、補間はおそらくかなり正確です)

ssizevec <- seq(10,250,by=5)
set.seed(101)
system.time(powvec <- sapply(ssizevec,powfun,m=5000))  ## 13 secs elapsed

これはかなり高速でありm=1e5、必要に応じて取得できる可能性がありますが、なぜそれほど正確な結果が必要なのかよくわかりません。電力曲線はかなり滑らかm=5000です...

長いシミュレーションを待ち焦がれている場合は、次のように置き換えることで、進行状況バーを印刷することもできsapply(ssizevec,powfun,m=5000)ますlibrary(plyr); aaply(ssizevec,.margins=1,powfun,.progress="text",m=5000)

最後に、カイ 2 乗値を直接サンプリングするか、分析力計算 (!) を行うことで、全体を大幅に高速化できると思います。これはループの最初の 2 行に相当すると思いますrchisq(m,df=ssize)*newvar/popvar。また、カイ 2 乗密度を直接数値計算することもできるかもしれません ...

system.time(powvec2 <- sapply(ssizevec,powfun,m=5000,sampchisq=TRUE))
## 0.24 seconds elapsed

（私はこれを試してみましたm=1e5.1から200までのサンプルサイズのすべての値でサンプリングしました... 24秒かかります...しかし、それでも不要かもしれないと思います。）

絵：

par(bty="l",las=1)
plot(ssizevec,powvec,type="l",xlab="sample size",ylab="power",
     xlim=c(0,250),ylim=c(0,1))
lines(ssizevec,powvec2,col="red")

ここに画像の説明を入力

于 2012-10-22T23:20:35.880 に答える
0

一般に、可読性/理解度ほど速度ではなく、可能な限りベクトル化を利用したいと考えています。

内部ループの内側に書き込むのはなぜですか(そして、計算も同様power[n]だと思います)。vecpals内側のループが実行された後、それは外側のループにあるべきではありませんか? 平方根の計算を両方のループの外に移動したい場合があります。

teststatvectとがpower行列 (明示的には 2 次元配列) として初期化され、ベクトル (または を使用して 1 次元配列として) として初期化されるのはなぜarrayですか? variance at 225前の行からのコメントの終わりだけですか? フォーマットを確認したい場合があります。（これは宿題ですか？）

ここでやろうとしているように見えることについてはreplicate、おそらくそれを呼び出す特定の関数を書くことによって、非常に便利な function を利用したいと思うかもしれません。

于 2012-10-22T22:53:36.147 に答える