記事の目的

分散も確率変数の場合の混合ガウスモデルの周辺化ギブスサンプリングをRで実装することです。
参考書籍: 佐藤一誠. ノンパラメトリックベイズ. 講談社. 2016
この書籍のアルゴリズム4.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)
# (1)
K <- 3
# (2)
z <- rep(1:K,N)[1:N]
# (3)
r <- 1
a0 <- 1
b0 <- 1

# (4)
m <- matrix(NA, nrow=K, ncol=D)
alpha <- 100
t <- matrix(NA, nrow=N, ncol=K)
n <- {}
max.iter <- 5
a <- 2*a0+(N-1)*D
x.av <- X %>% apply(2, mean)
x.av <- matrix(rep(x.av,each=N-1), nrow=N-1, ncol=D)
x.k.i <- matrix(NA, nrow=K, ncol=D)

for(s in 1:max.iter){
  # (ⅰ)
  for(i in 1:N){
    X.i <- X[-i,]
    rownames(X.i) <- 1:(N-1)
    z.i <- z[-i]
    x.av.i <- X.i %>% apply(2, mean)
    for(k in 1:K){
      # x.k\i
      x.k.i[k,] <- X.i[z.i==k,] %>%  apply(2, mean)
      # n\i
      n[k] <- z.i[z.i==k] %>% length()
      # m\i
      m[k,] <- n[k]*x.k.i[k,]/(n[k]+r)
    }
    #b\i
    b <- 2*b0+(unlist(X.i-x.av) %*% unlist(X.i-x.av))+(N-1)*r/N*(unlist(x.av.i) %*% unlist(x.av.i))
    for(k in 1:K){
      t[i,k] <- dmvt(as.numeric(X[i,]), delta=m[k,], df=a, sigma=(1+1/N)*b[1]*diag(D))*((n[k]+alpha)/(sum(n)+alpha*K))
    }
    z_prob <- t[i,]/ sum(t[i,])
    z[i] <- rmultinom(1, 3, z_prob) %>% which.max()
  }
}

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)

補足

コードの雰囲気は良い感じだったんですが、全然クラスタリングできてない…
おそらくt分布の式が違うんですが、対処法がわかりしだい修正します、、、
誰か間違いが分かる人いたら教えて欲しいです😥