Rコードの最適化例：混合正規乱数の発生コード

Rコードの最適化例：混合正規乱数の発生コード (Rコード最適化のコツと実例集 に戻る)

論理ベクトルの使用法の見本

舟尾さんによる正規分布 N(-1, 1), N(2,1) を確率 0.6 で混合した混合正規分布確率変数の発生コード

出典は Smoothing Techniques with Implementation in S. （Haerdle, W. (1991). New York: Springer-Verlag.）です．

generator <- function(n) {
   data1 <- rnorm(n)-1           # N(-1,1) に従う乱数
   data2 <- rnorm(n)+2           # N(2, 1) に従う乱数
   data3 <- (runif(n) <= 0.6)    # 確率 0.6 で 1 となる論理ベクトル
   data1*data3 + data2*(1-data3) # 確率 0.6 で N(-1,1), 確率0.4で N(2,1) 乱数を採択
}

なんでも掲示版に投稿された改良版(2n 個の rnorm 乱数のうち実際使うのは n 個だけという無駄を省く)

## -1 + 3*(runif(n) > 0.6) は確率 0.6 で 2、確率 0.4 で -1 となるベクトル
generator2 <- function(n) { 
  (rnorm(n)-1)+3*(runif(n) > 0.6) 
}

応用問題：二つの任意正規分布 N(mu1, sd1^2), N(mu2, sd2^2) を確率 p1:(1-p1) で混合

# n: 発生する混合正規乱数の総数
# p1: 最初の正規分布の比率
# N1: 最初の正規分布の平均、標準偏差のベクトル c(mu1, sd1)
# N2: 2番目の正規分布の平均、標準偏差のベクトル c(mu2, sd2)
rmixnorm <- function(n, p1, N1, N2) {
   P <- (runif(n) < p1)
   (P*N1[2] + (1-P)*N2[2])*rnorm(n) + (P*N1[1] + (1-P)*N2[1])
}

応用問題2：見掛けを少しすっきりさせる

rmixnorm2 <- function(n, p1, N1, N2) {
   P <- (runif(n) < p1)
   NN <-  N1%o%P + N2%o%(1-P) # 外積演算を使い改良(悪?)
   NN[2,]*rnorm(n) + NN[1,]
}

何でも掲示版記事より

結局 N1 に従う乱数の個数 n1 が決まれば良いと喝破(なるほど)

rmixnorm3 <- function(n, p1, N1, N2) {
   n1 <- sum(runif(n) < p1)
   c( N1[2]*rnorm(n1) + N1[1], N2[2]*rnorm(n-n1)+N2[1] )
}

気持ちが悪ければ結果を sample 関数で並べ変える。

更に rnorm の平均、標準偏差引数を使えば済むとの指摘(然り)

rmixnorm4 <- function(n, p1, N1, N2) {
   n1 <- sum(runif(n) < p1)
   c( rnorm(n1, m=N1[1], s=N1[2]), rnorm(n-n1, m=N2[1], s=N2[2]) )
}

実行速度比較(100 回の平均)

>  mean(sapply(1:100, function(i) system.time(generator(1e6))[1]))
[1] 2.7396
> mean(sapply(1:100, function(i) system.time(generator2(1e6))[1])) 
[1] 1.5031

> mean(sapply(1:100, function(i) system.time(rmixnorm(1e6, 0.4, c(-1,1), c(3,2)))[1]))
[1] 2.3924
>  mean(sapply(1:100, function(i) system.time(rmixnorm2(1e6, 0.4, c(-1,1), c(3,2)))[1])) 
[1] 2.8277
> mean(sapply(1:100, function(i) system.time(rmixnorm3(1e6, 0.4, c(-1,1), c(3,2)))[1]))
[1] 1.4865
> mean(sapply(1:100, function(i) system.time(rmixnorm4(1e6, 0.4, c(-1,1), c(3,2)))[1])) 
[1] 1.4215

演習問題

任意個数の正規分布を任意比率で混合した混合正規分布に従う乱数を発生させるうまいコードを書け。ヒント：rmixnorm4 が参考になる。

最新の20件