廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

混合正規乱数の生成

混合正規乱数を生成する R のコードを書け

問:

任意個数の正規分布を任意比率で混合した混合正規分布に従う乱数を発生させるうまいコードを書け。
Rコードの最適化例:混合正規乱数の発生コード - RjpWiki より)

答:

rmixnorm <- function(n,mu, sigma,weight){
  stopifnot( length(mu) == length(sigma) & length(sigma)==length(weight))
  Y <-sample(length(weight), n, replace=TRUE, prob=weight)
  (X <-rnorm(n, mu[Y], sigma[Y]))
}

実行例:

mu <- c(1, 2, 5, 4)
sigma <- sqrt(c(0.01, 0.5, 0.02, 0.01))
weight <- c(0.1, 0.6, 0.2, 0.1)
X <- rmixnorm(10000,mu,sigma,weight)

hist(X, breaks="scott", prob=TRUE,main=NULL,col="gray80", border="gray20")
curve(weight[1]* dnorm(x,mu[1],sigma[1])+
  weight[2] *dnorm(x,mu[2], sigma[2])+
  weight[3] *dnorm(x,mu[3],sigma[3])+
  weight[4] *dnorm(x,mu[4],sigma[4]),
  col="red4", lwd=2,
  min(X), max(X), n=1000, add=TRUE)
#n: the number of x values at which to evaluate.
#デフォルトでは n = 101 でジャギーなカーブになる

f:id:abrahamcow:20140621032511p:plain

解説:

An Introduction to Statistical Computing: A Simulation-based Approach (Wiley Series in Computational Statistics)

An Introduction to Statistical Computing: A Simulation-based Approach (Wiley Series in Computational Statistics)

↑この本にのっていた sample 関数の使い方がうまいなーと思ったので、まんま転載しました。

混合分布(mixture distribution)F_\theta とは、確率分布関数  F_1, \ldots, F_k、重み  \theta_1, \ldots, \theta_k(ただし \sum^k_{i=1} \theta _i =1) に対して、
 \displaystyle F_\theta = \sum^{k}_{i=1}\theta_i F_i(x)
で与えられる分布です。