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

No.21761　多項分布に従う乱数発生について　　【看護研究屋】　2015/09/04(Fri) 10:33

いつも勉強させていただいております。

シミュレーションを行うために，n=100万ほどのデータを発生させたいと考えています。
シミュレーションの設定は，次のとおりです。
・個人毎に曝露される強度が異なる（0～1の一様分布），2種類の曝露がある
・アウトカム（その個人がどうなるか）が3種類あり（1or2or3），曝露によってリスクが変わる
　・z1への曝露が多いと，1にはなりにくくなる
　・z2への曝露が多いと，3になりやすくなる

そこで，次のようにプログラムを書いたのですが，n=1万ならともかく，100万とするとfor文のところで5時間ほどかかってしまい，どうにか短縮できないかご相談する次第です。

n <- 10000
d <- data.frame(ID=matrix(rep(1:n,1),nrow=n))
d$x1 <- c(runif(n, min=0, max=1))
d$x2 <- c(runif(n, min=0, max=1))
d$p1 <- 0.5 - d$x1 / 10 - d$x2 / 10
d$p2 <- 0.4 + d$x1 / 20 - d$x2 / 10
d$p3 <- 0.1 + d$x1 / 20 + d$x2 / 5

for(i in 1:n){
temp <- rmultinom(1, size=1, prob=c(d$p1[i], d$p2[i], d$p3[i]))
d$y1[i] <- temp[1]
d$y2[i] <- temp[2]
d$y3[i] <- temp[3]
}

d$y <- as.factor(d$y1 + 2 * d$y2 + 3 * d$y3)
summary(d)

行ごとに処理するのではなく，一度に処理できないか検討した（調べた）のですが，調べきれませんでした。ご指導何卒宜しくお願い申し上げます。

No.21762　Re: 多項分布に従う乱数発生について　　【青木繁伸】　2015/09/04(Fri) 14:27

あなたのプログラムも，以下のようにすると3，4倍速くなりますね。
n=1000000 でも 50 秒で計算できます。このようにしたプログラムの計算時間は，n に比例します。
あなたのプログラムはその都度使い捨てられる中間変数 temp を使っているので，頻繁にガーベッジコレクションが発生しているものと見られます。このため計算時間は n に比例する以上に必要になっている（処理時間のほとんどがガーベッジコレクションになっているのでしょう。教訓：無用な中間変数は使わない。ループの中で余計な処理をしない。）。
n <- 1e+06
d <- data.frame(ID = matrix(rep(1:n, 1), nrow = n))
d$x1 <- c(runif(n, min = 0, max = 1))
d$x2 <- c(runif(n, min = 0, max = 1))
d$p1 <- 0.5 - d$x1/10 - d$x2/10
d$p2 <- 0.4 + d$x1/20 - d$x2/10
d$p3 <- 0.1 + d$x1/20 + d$x2/5

d2 <- matrix(0, n, 3)
for (i in 1:n) {
	d2[i, ] <- rmultinom(1, size = 1, prob = c(d$p1[i], d$p2[i], d$p3[i]))
}
d$y1 <- d2[, 1]
d$y2 <- d2[, 2]
d$y3 <- d2[, 3]

d$y <- as.factor(d$y1 + 2 * d$y2 + 3 * d$y3)
summary(d)

No.21763　Re: 多項分布に従う乱数発生について　　【看護研究屋】　2015/09/04(Fri) 22:17

ご教示ありがとうございます！
おかげさまで，速くすることができました。
ガーベッジコレクションの用語・概念を初めて知りました。
教訓もいただき，今後にも活かしたいと思います。
ほんとうにありがとうございました。