4

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

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

シミュレーションの高速化に関するコメントをいただければ幸いです。

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))
    }
}

データセット内の 7135 個の観測値のそれぞれについて、100 個のシミュレートされた栄養値が作成され、元の測定レベルに逆変換されます (シミュレーションは、BoxCox 変換された栄養値に対する非線形混合効果モデルの結果を使用しています)。

forループは非効率的であると読んだので、ループを使用しないことをお勧めしますが、それらを代替として使用するためRのオプションについて十分に理解していません。スタンドアロン マシンで実行されている場合、コードの変更方法に関する推奨事項に影響する場合、通常、これは Windows 7 バリアントを実行する標準の Dell タイプのデスクトップになります。applyR

更新: テストのためにこれを再現するには、 Lambda.Value=0.4 およびMale.Resid.Var=12.1029420429778 でMale.Distrib$stddev_u2あり、すべての観測値で一定の値です。

str(Male.Distrib)は

'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

NaN更新 2:結果を引き起こしている関数の行は

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

支援とコメント、そして応答の速さに感謝します。

更新: @Ben Bolker はtemp、NaN の問題を引き起こしているのは負の値であることは正しいです。いくつかのテストでこれを見逃しました（値のみが返されるように関数をコメントアウトしtemp、結果の data frame を呼び出した後Test）。このコードはNaN問題を再現します:

> min(Test)
[1] -2.103819
> min(Test)^(1/Lambda.Value)
[1] 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ステートメントは、実行する代わりのものです。
4

1 に答える 1

4

2 つの最大の速度の問題に対処するアプローチを次に示します。

  1. 観察 ( i) をループする代わりに、それらをすべて一度に計算します。
  2. jMC レプリケーション ( )をループする代わりに、この目的のため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 の追加などを行うことができますが、これは主要な計算部分のアイデアです。幸運を！

于 2012-01-25T20:25:24.840 に答える