Rで混合正規分布を描く。

Rで混合正規分布をプロットしようとして、意外に時間がかかったので、ソースコードをメモ。

ちなみに多変量正規分布を扱うパッケージは、mvtnorm
というのがあったのでこいつを使用した。

library(mvtnorm)

#MixtureModelDens：混合分布を描く元になる密度関数を描く
#from,toの範囲内で、length数で、混合分布を描く
#m1,m2：2つの混合分布の平均となる2次元ベクトル
#ｓ1,ｓ2：2つの混合分布の標準偏差となる2次元行列
#pi1,pi2：2つの混合分布の選択率
MixtureModelDens<-function(from, to, length, m1, m2, s1, s2, pi1, pi2){

 #計算の範囲を算出する。
　x1<-seq(from=from, to=to, length=length)
 x2<-x1

 #密度計算関数の作成する。
 f1<-function(x1, x2){
	dmvnorm(matrix(c(x1, x2), ncol=2), mean=m1, sigma=sigma1)
 }
 f2<-function(x1, x2){
 	dmvnorm(matrix(c(x1, x2), ncol=2), mean=m2, sigma=sigma2)
 }
　#個々の正規分布を計算する。
 z1<-outer(x1, x2, f1)
 z2<-outer(x1, x2, f2)
 
 pi1*z1+pi2*z2
}

s1<-matrix(c(1,0.2,0.2,1), ncol=2)
s2<-matrix(c(1,0.4,0.4,1), ncol=2)

x<-seq(from=-7, to=7, length=140)
y<-seq(from=-7, to=7, length=140)
z<-MixtureModelDens(-7, 7, 140, c(1,1), c(-4, -2), s1, s2, 0.4, 0.6)

op<-par(bg = "white")
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")

結果は以下の画像

ちなみに参考にしたのは以下のURL
http://cse.niaes.affrc.go.jp/minaka/R/R-binormal.html