Rでシミュレーションを行おうと思って、「群変数が共通で、データを乱数として r セット生成する」スクリプトをループを使わずに書こうと思ったら、なかなかうまくいかない。
↑を逆(「データが共通で、群変数を乱数として r セット生成する」)にするとリサンプリングやブートストラップも出来そう。
set.seed(090930) r <- 2 m <- c(2, 3, 4) s <- c(.1, .1, .1) n <- c(2, 3, 4) d <- do.call("rbind", mapply(function(n, m, s) matrix(rnorm(r * n, mean = m, sd = s), nrow=n, ncol=r), n, m, s, SIMPLIFY = FALSE) ) g <- rep(1:length(n), n) d; g apply(d, 2, function(d, g) tapply(d, g, mean), g)
> d; g [,1] [,2] [1,] 2.003117 2.179617 [2,] 1.938691 1.998216 [3,] 2.851668 3.109537 [4,] 2.982875 3.088262 [5,] 3.051431 3.047164 [6,] 3.949051 4.144669 [7,] 4.071078 3.878844 [8,] 3.892767 4.207359 [9,] 3.963631 3.989085 [1] 1 1 2 2 2 3 3 3 3 > apply(d, 2, function(d, g) tapply(d, g, mean), g) [,1] [,2] 1 1.970904 2.088917 2 2.961991 3.081655 3 3.969132 4.054989
Rのデータハンドリングの難しさが身にしみました。
以下のシミュレーションに使う計算もなるべくapply系関数を駆使して書いておきます。
snow packageのparApply系関数に変えると最近のマルチコアCPUだと計算時間が 1/石の数 ぐらいになるはずです。
パッケージ開発者は、parMapply関数がそのうち実装ってだいぶ前に言ってるけどどうなんだろう。
追記
09/10/21
mapply の部分にオプション "SIMPLIFY = FALSE" を追加 (修正)。
等標本数の場合これを入れないとエラーが起きてしまう。
