単変量混合正規分布モデルをデータにフィットさせる

外国語教育研究では2つの山とか3つの山とかの分布になっているのを見ることがある。こういうときは，混合分布モデルをデータにフィットさせるといいかもだ。Rではmclustもいいけど，mixtoolsというパッケージがある。

#準備
library(mixtools)

#数値例の作成
#ここでは混合比（λ）が1:2，母平均（μ）が50, 120，母標準偏差（σ）が10, 20が正解
set.seed(92)
huta<-c(rnorm(100,50,10),rnorm(200,120,20))

#可視化して確認
par(mfrow=c(1,2))
hist(huta,col="lightgray",xlab="Value",main="")
plot(ecdf(huta))

#mixtoolsのEMアルゴリズム関数を使って母数を推定
fit<-normalmixEM(huta)
summary(fit)

summary of normalmixEM object:
           comp 1    comp 2
lambda   0.669647  0.330353
mu     117.657305 50.067067
sigma   20.090447  9.126634
loglik at estimate:  -1428.916 

#ここではだいたい正解だった
#対数尤度
fit$loglik
#AIC　（自由母数は５個でいいよね？）
fitaic<--2*fit$loglik+2*5
fitaic

#誤差をブートストラップで
boot.se(fit)

#ヒストグラムに確率密度関数を合わせる
hist(huta,col="lightgray",xlab="Value",main="",freq=F,ylim=c(0,.015),breaks=20)

#その前に推定した母数をもつ混合正規分布の確率密度関数を定義する
dnormmix<-function(x){
y<-fit$lambda[1]*dnorm(x,fit$mu[1],fit$sigma[1])+fit$lambda[2]*dnorm(x,fit$mu[2],fit$sigma[2])
y
}

x<-seq(0,200,0.1)
lines(x,dnormmix(x),lwd=2,col=2)

うむ。間便な混合正規分布の確率密度関数も書いておく。これで好きにお絵かきできる。

dnormmix<-function(x,lambda, mu1, mu2, sigma1, sigma2){
y<-lambda*dnorm(x,mu1,sigma1)+(1-lambda)*dnorm(x,mu2,sigma2)
y
}

（なぜか混合分布モデルは自分の青春だった，って感じがする）