記事の目的

平均も分散も確率変数の場合の混合ガウスモデルのギブスサンプリングをRで実装することです。
参考書籍: 佐藤一誠. ノンパラメトリックベイズ. 講談社. 2016
この書籍のアルゴリズム4.3に対応します
※目次2, 4のコードは前回と同様です。

目次

1. モデルの説明
 2. 使用データ
 3. 分散固定の混合ガウスモデルのギブスサンプリング
 4. クラスタリング確認

1. モデルの説明

2. 使用データ

> X <- iris[,1:4]
> D <- ncol(X)
> N <- nrow(X)
> head(X)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
>

3. 分散固定の混合ガウスモデルのギブスサンプリング

library(dplyr)
library(mvtnorm)
library(MCMCpack)
# (1)
K <- 3
# (2)
u <- matrix(rep(apply(X, 2, mean),each=K), nrow=K, ncol=D)
# (3)
r <- 1
a0 <- 1
b0 <- 1
# (4)
pi <- rep(1/K, K)
t <- 1

# (5)
alpha <- 1/K
max.iter <- 30
Norm <- {}
z <- {}
n <- {}
x.k <- matrix(NA, nrow=K, ncol=D)
tmp.sum2 <- {}
a <- a0+N*D/2

set.seed(100)
for(s in 1:max.iter){
  tmp.sum1 <- 0
  # (ⅰ)
  for(i in 1:N){
    for(k in 1:K){
      Norm[k] <- dmvnorm(X[i,], u[k,], diag(D)/sqrt(t))*pi[k]
    }
    z_prob <- Norm/sum(Norm)
    z[i] <- rmultinom(1, 3, z_prob) %>% which.max()
  }
  # (ⅱ)
  for(k in 1:K){
    n[k] <- z[z==k] %>% length()
    x.k[k,] <- X[z==k,] %>%  apply(2, mean)
    u[k,] <- rmvnorm(1, n[k]*x.k[k,]/(n[k]+r), diag(D)/sqrt(t*(n[k]+r)))
  }
  # (ⅲ)
  for(k in 1:K){
    for(i in 1:N){
      if(z[i]==k){
        tmp.sum1 <- tmp.sum1 + unlist(X[i,]-x.k[k,]) %*% unlist(X[i,]-x.k[k,])
      }
    }
    tmp.sum2[k] <- tmp.sum1/2+n[k]*r/(2*(r+n[k]))*(t(x.k[k,]) %*% (x.k[k,]))
  }
  b <- b0 + sum(tmp.sum2)
  t <- rgamma(1,a,b)
  # (ⅳ)
  alpha <- alpha + n
  pi <- rdirichlet(1, alpha)
}

4. クラスタリング確認

library(cluster)
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)