Rで分散固定の場合の混合ガウスモデルの実装 (ギブスサンプリング)

記事の目的

分散固定の場合の混合ガウスモデルをRで実装します。
for文を極力使用せず、早く学習するようにしています。
参考: ノンパラメトリックベイズ 点過程と統計的機械学習の数理

目次

No.	目次
1	モデルの説明
2	データとライブラリ
3	実装
4	確認

1. モデルの説明

2. データとライブラリ

データirisを使います。

X <- iris[,1:4]
D <- ncol(X)
N <- nrow(X)
library(mvtnorm)
library(cluster)

3. 実装

#(1)Kを求める
K <- 3
#(2)muを乱数で初期化
set.seed(100)
mu <- matrix(rep(apply(X,2,mean),each=K)+rnorm(K*D,0,0.2), nrow=K)
#(3)(ⅰ)(ⅱ)を繰り返す
max.iter <- 30
for(s in 1:max.iter){
  #(ⅰ)zのサンプリング
  tmp <- apply(mu, 1, function(x) dmvnorm(X, x, diag(D)))
  z <- apply(tmp, 1, function(x) which.max(rmultinom(1,3,x)))
  #(ⅱ)muのサンプリング
  n <- tapply(z, z, length)
  x.k <-  apply(X, 2, function(x) tapply(x, z, mean))
  mu <- t(apply(cbind(x.k, n), 1,
        function(x) rmvnorm(1, x[D+1]*x[1:(D)]/(x[D+1]+1), diag(D)/(x[D+1]+1))))
}

4. 確認

左が正解で、右が実装の結果です。

par(mfrow=c(1,2))
clusplot(X, iris[,5], color=TRUE, shade=FALSE, labels=4, lines=0)
clusplot(X, z, color=TRUE, shade=FALSE, labels=4, lines=0)