読者です 読者をやめる 読者になる 読者になる

箱形に乱数を生成する

R

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" を追加 (修正)。
等標本数の場合これを入れないとエラーが起きてしまう。