Rで無限混合ガウスモデルへの道.4(未完成)

記事の目的

無限混合ガウスモデルの周辺化ギブスサンプリングをRで実装することです。
参考書籍: 佐藤一誠. ノンパラメトリックベイズ. 講談社. 2016

目次

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 <- 1
# (2)
z <- rep(1:K,N)[1:N]
# (3)
m0 <- rep(0, D) 
r <- 1
a0 <- 1
b0 <- 1

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

set.seed(100)
for(s in 1:max.iter){
  t <- matrix(NA, nrow=N, ncol=K)
  m <- matrix(NA, nrow=K, ncol=D)
  x.av <- matrix(rep(x.av,each=N-1), nrow=N-1, ncol=D)
  x.k.i <- matrix(NA, nrow=K, ncol=D)
  # (ⅰ)
  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])/(N-1+alpha))
    }
    t1[i] <- dmvt(as.numeric(X[i,]), delta=m0, df=2*a0, sigma=(1+1/r)*b0*diag(D))*(alpha/(N-1+alpha))
    z_prob <- append(t[i,], t1[i])/ sum(append(t[i,], t1[i]))
    z[i] <- rmultinom(1, 1, z_prob) %>% which.max()
  }
  K <- z %>% unique() %>% length()
}

4. クラスタリング確認

前回同様t分布が違うと思います。
今回はエラーで終わりました間違いが分かり次第修正します、、、