より効率的なモンテカルロシミュレーションループの作成方法R

Question

この練習の目的は、栄養摂取量の母集団分布を作成することです。以前のデータには繰り返し測定値がありましたが、これらは削除されているため、各行はデータフレーム内の一意の人物です。

私はこのコードを持っています。このコードは、少数のデータフレーム行でテストしたときに非常にうまく動作します。 7135行すべてに対して、非常に遅いです。時間を計ろうとしましたが、マシンの経過時間が15時間のときにクラッシュしました。 system.timeの結果はTiming stopped at: 55625.08 2985.39 58673.87であった。

私はシミュレーションの高速化上の任意のコメントをいただければ幸いです：私のデータセットで7135個の観測値のそれぞれについて、

Male.MC <-c() 
for (j in 1:100)   { 
for (i in 1:nrow(Male.Distrib)) { 
    u2  <- Male.Distrib$stddev_u2[i] * rnorm(1, mean = 0, sd = 1) 
    mc_bca <- Male.Distrib$FixedEff[i] + u2 
    temp  <- Lambda.Value*mc_bca+1 
    ginv_a <- temp^(1/Lambda.Value) 
    d2ginv_a <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2)) 
    mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var/2 
z <- data.frame(
    RespondentID = Male.Distrib$RespondentID[i], 
    Subgroup  = Male.Distrib$Subgroup[i], 
    mc_amount = mc_amount, 
    IndvWeight = Male.Distrib$INDWTS[i]/100 
    ) 

Male.MC <- as.data.frame(rbind(Male.MC,z)) 
    } 
} 

を、100のシミュレートされた栄養素の値が作成され、その後、元の測定レベル（シミュレーションに変換BoxCox変換栄養素値に対する非線形混合効果モデルの結果を使用している）。

私は、彼らがRに非効率的であるが、私は、代替として、それらを使用するようにapplyに基づくオプションについて十分に理解していないことを読んで私は、forループを使用しないでしょう。 Rはスタンドアローンのマシンで実行されていますが、通常、これはWindows 7の亜種を実行する標準のDellタイプのデスクトップで、コードを変更する方法の推奨事項に影響する場合に使用します。

更新：テストのためにこれを再現するために、 Lambda.Value = 0.4 Male.Resid.Var = 12.1029420429778とMale.Distrib$stddev_u2はすべての観測にわたって一定の値です。

​​は

'data.frame': 7135 obs. of 14 variables: 
$ RndmEff  : num 1.34 -5.86 -3.65 2.7 3.53 ... 
$ RespondentID: num 9966 9967 9970 9972 9974 ... 
$ Subgroup : Ord.factor w/ 6 levels "3"<"4"<"5"<"6"<..: 4 3 2 4 1 4 2 5 1 2 ... 
$ RespondentID: int 9966 9967 9970 9972 9974 9976 9978 9979 9982 9993 ... 
$ Replicates : num 41067 2322 17434 21723 375 ... 
$ IntakeAmt : num 33.45 2.53 9.58 43.34 55.66 ... 
$ RACE  : int 2 3 2 2 3 2 2 2 2 1 ... 
$ INDWTS  : num 41067 2322 17434 21723 375 ... 
$ TOTWTS  : num 1.21e+08 1.21e+08 1.21e+08 1.21e+08 1.21e+08 ... 
$ GRPWTS  : num 41657878 22715139 10520535 41657878 10791729 ... 
$ NUMSUBJECTS : int 1466 1100 1424 1466 1061 1466 1424 1252 1061 1424 ... 
$ TOTSUBJECTS : int 7135 7135 7135 7135 7135 7135 7135 7135 7135 7135 ... 
$ FixedEff : num 6.09 6.76 7.08 6.09 6.18 ... 
$ stddev_u2 : num 2.65 2.65 2.65 2.65 2.65 ... 

head(Male.Distrib)ある

RndmEff RespondentID Subgroup RespondentID Replicates IntakeAmt RACE INDWTS TOTWTS GRPWTS NUMSUBJECTS TOTSUBJECTS FixedEff stddev_u2 
1 1.343753   9966  6   9966  41067 33.449808 2 41067 120622201 41657878  1466  7135 6.089918 2.645938 
2 -5.856516   9967  5   9967  2322 2.533528 3 2322 120622201 22715139  1100  7135 6.755664 2.645938 
3 -3.648339   9970  4   9970  17434 9.575439 2 17434 120622201 10520535  1424  7135 7.079757 2.645938 
4 2.697533   9972  6   9972  21723 43.340180 2 21723 120622201 41657878  1466  7135 6.089918 2.645938 
5 3.531878   9974  3   9974  375 55.660607 3 375 120622201 10791729  1061  7135 6.176319 2.645938 
6 6.627767   9976  6   9976  48889 91.480049 2 48889 120622201 41657878  1466  7135 6.089918 2.645938 

アップデート2：NaN結果を引き起こしている機能の行は、彼らの支援のためのみんなに

d2ginv_a <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2)) 

おかげでコメント、およびspe応答のed。

更新：@Ben Bolkerは、NaNの問題の原因となる値が負の値tempであるという点には間違いありません。私はいくつかのテスト（関数をコメントアウトしてtempの値だけが返され、結果データフレームTestが呼び出されたのでこれを逃しました）。

> min(Test) 
[1] -2.103819 
> min(Test)^(1/Lambda.Value) 
[1] NaN 

しかし、値としての値を入れ、その後、同じことを実行している計算は私に結果を与えるので、手計算を行うとき、私はこれを逃した（？）：このコードはNaN問題を再現

> -2.103819^(1/Lambda.Value) 
[1] -6.419792 

私はベクトル化を使用している（私が思う）働くコードを持っており、それは驚くほど高速です。誰かがこの問題を抱えている場合に備えて、私は以下の作業コードを掲示しています。私は< 0計算の問題を防ぐために最小値を加えなければならなかった。助けてくれた皆様、コーヒーに感謝します。私はrnormの結果をデータフレームに入れようとしましたが、実際には遅くなりました。このようにして、cbindを使うのは本当に速いです。Male.Distribは私の7135観測データの完全なデータフレームですが、このコードは私が前に投稿したカットダウン版（テストされていない）で動作するはずです。

Min_bca <- ((.5*min(Male.AddSugar$IntakeAmt))^Lambda.Value-1)/Lambda.Value 
Test <- Male.Distrib[rep(seq.int(1,nrow(Male.Distrib)), 100), 1:ncol(Male.Distrib)] 
RnormOutput <- rnorm(nrow(Test),0,1) 
Male.Final <- cbind(Test,RnormOutput) 
Male.Final$mc_bca <- Male.Final$FixedEff + (Male.Final$stddev_u2 *  Male.Final$RnormOutput) 
Male.Final$temp  <- ifelse(Lambda.Value*Male.Final$mc_bca+1 > Lambda.Value*Min_bca+1, 
          Lambda.Value*Male.Final$mc_bca+1, Lambda.Value*Min_bca+1) 
Male.Final$ginv_a <- Male.Final$temp^(1/Lambda.Value) 
Male.Final$d2ginv_a <- ifelse(0 > (1-Lambda.Value)*Male.Final$temp^(1/Lambda.Value-2), 
          0, (1-Lambda.Value)*Male.Final$temp^(1/Lambda.Value-2)) 
Male.Final$mc_amount <- Male.Final$ginv_a + Male.Final$d2ginv_a * Male.Resid.Var/2 

日のレッスン：あなたは私はあなたがmax()使用することはできません以前

をしようとしていた何をしようとした場合

• 分布関数は、ループ内でリサンプリングされていないよう私は2つの値から最大値を求めていましたが、カラムから最大値を返していました。 ifelseステートメントは、これを置き換えるものです。ここで

Answer 1

2つの最大のスピードの問題に対処するアプローチです：

1. 代わり観測（i）をループで、我々は一度にすべてを計算します。
2. MCレプリケーション（j）をループする代わりに、replicateを使用します。これは、この目的のための簡略化されたapplyです。

まず、データセットを読み込み、実行中の機能を定義します。

Male.Distrib = read.table('MaleDistrib.txt', check.names=F) 

getMC <- function(df, Lambda.Value=0.4, Male.Resid.Var=12.1029420429778) { 
    u2  <- df$stddev_u2 * rnorm(nrow(df), mean = 0, sd = 1) 
    mc_bca <- df$FixedEff + u2 
    temp  <- Lambda.Value*mc_bca+1 
    ginv_a <- temp^(1/Lambda.Value) 
    d2ginv_a <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2)) 
    mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var/2 
    mc_amount 
} 

次に、私たちはそれを複数回複製します。

> replicate(10, getMC(Male.Distrib)) 
     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] 
[1,] 36.72374 44.491777 55.19637 23.53442 23.260609 49.56022 31.90657 25.26383 25.31197 20.58857 
[2,] 29.56115 18.593496 57.84550 22.01581 22.906528 22.15470 29.38923 51.38825 13.45865 21.47531 
[3,] 61.27075 10.140378 75.64172 28.10286 9.652907 49.25729 23.82104 31.77349 16.24840 78.02267 
[4,] 49.42798 22.326136 33.87446 14.00084 25.107143 25.75241 30.20490 33.14770 62.86563 27.33652 
[5,] 53.45546 9.673162 22.66676 38.76392 30.786100 23.42267 28.40211 35.95015 43.75506 58.83676 
[6,] 34.72440 23.786004 63.57919 8.08238 12.636745 34.11844 14.88339 21.93766 44.53451 51.12331 

次に、再フォーマット、IDの追加などができますが、これは主な計算部分の考えです。がんばろう！