「統計学関連なんでもあり」の過去ログ--- 043

No.12927　モンテカルロシミュレーションのグラフ　　【キムキム】　2010/06/26(Sat) 18:37

初めて投稿させていただきます。1週間以上悩んでいますがなかなか前に進めずこちらを見つける次第となりました。統計学というよりはRの使い方，になっておりますが私にはここ以外に御相談できる場所がなく，どうかお答え頂けますと幸いです。

Rを用いたモンテカルロシミュレーションのために薬物血中濃度の曲線を描いています。
薬物の消失定数KeとクリアランスCLに乱数を与えることはできたのですが，複数（例えば20ぐらい）のグラフをいっぺんに表示したり，いっぺんにAUCを計算したりする方法がわかりません。for を用いた繰り返しでグラフ表示ができるのかなとやってみた例が以下ですがうまくいきませんでした。

Ke とCLに乱数を与えたデータフレームをzとしました。
ちなみにDoseは1000(mg) ，infT は投与時間で1.0 時間という意味です。

> Y <- function(x) Dose/CL/infT*(1-exp(-Ke*infT))*exp(-Ke*x)

> Dose <- 1000 
> infT <- 1.0 

> z
   ID Clearance        KE
1   1  71.72329 0.4405448
2   2  80.25860 0.4263393
3   3  75.42981 0.4274200
4   4  46.40765 0.4065925
5   5  63.61894 0.3784054
6   6  56.81234 0.4133871
7   7  69.48724 0.3908536
8   8  74.37029 0.4044656
9   9  74.31195 0.4312741
10 10  30.63537 0.3946689
11 11  73.83538 0.4238592
12 12  67.21356 0.4261859
13 13  54.08096 0.4420179
14 14  69.37479 0.4334032
15 15  58.20908 0.4176028
16 16  44.66195 0.4247530
17 17  57.42262 0.4184082
18 18  42.72584 0.4140849
19 19  87.49260 0.3959427
20 20  52.23639 0.4207313

> myloop
function() {
    a <- 0
    for (i in 1:20) {
        a <- a+1
        CL <- z[a, 2]
        Ke <- z[a, 3]
        curve(Y, col=4, add=TRUE)
    }
}
> myloop()

でグラフが一本しか表示されません。curve までの内容がループされていないのでしょうか？素人の質問で誠に申し訳ありませんが3日以上悩んでこの結果までしか到達できませんでした。
どうぞご教授よろしくお願い申し上げます。

No.12928　Re: モンテカルロシミュレーションのグラフ　　【青木繁伸】　2010/06/26(Sat) 20:19

R コンソールにコピー・ペーストして実行できるものを提示してください。

あれこれ，文句を言われて，実行することができません。
　　オブジェクト 'CL' がありません
　　オブジェクト 'Ke' がありません
　　まだ plot.new が呼ばれていません

適当に直して動いたとしても，あなたがやっていることを再現できていないと思うので，そこからスタートしてああだこうだアドバイスしても無意味ですよね。

グラフが一本しか描けないというのは，R でグラフを描画するときの基本概念が把握できていないためでしょうね。

curv は，基本的には plot を使っていて，plot が実行されるたびに，「新しいグラフィック」が描かれます。重ね描きするには，二回目以降の描画は，lines とか plots とかを使うとか plot.new とか par(new=) を使う。

例えば以下のように？
Y は myloop の中で定義する
図を描く前に par(new=TRUE) を実行する
そのようにすれば一応重ね描きができるようだけど，それがあなたの望んでいる正しいグラフなのかどうかは，私にはわからない。
Dose <- 1000 
infT <- 1.0
myloop <- function() {
    Y <- function(x) Dose/CL/infT*(1-exp(-Ke*infT))*exp(-Ke*x)
    old <- par(new=TRUE)
    a <- 0
    for (i in 1:20) {
        a <- a+1
        CL <- z[a, 2]
        Ke <- z[a, 3]
        curve(Y, col=4, add=TRUE)
    }
    par(old)
}
myloop()

No.12934　Re: モンテカルロシミュレーションのグラフ　　【キムキム】　2010/06/27(Sun) 02:27

早々のご返信誠にありがとうございました。
記載からもう少し自分なりに勉強させていただきます。
私にとりましては先生の記載が求めていた情報でございました。

どうもありがとうございました。