廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

WBICで混合多項分布のクラスタ数を決めてみる(ギブスサンプリング)

今日の川柳

Stan を使って WBIC を計算する例は、
WAICとWBICを事後分布から計算する - StatModeling Memorandum

Stan を使えばWBICの計算は簡単なのですが、場合によってはギブスサンプリングなどのアルゴリズムを自分で導出したいこともあると思います。

WBIC の計算には逆温度 \beta の事後分布に従う乱数を作る必要があります。

これを作るには、要は尤度関数を \beta 乗すればいいのだと理解しています(まちがってたらすみません。)

今回は混合多項分布モデルを考えます。

各サンプルがどのクラスタから得られたかを示すインジケータ変数 z_{nk} が得られているとすると、混合多項分布の完全な尤度関数は

 \prod_{n=1}^{N} \prod_{l=1}^{L} \phi_l^{z_{nl}} p_{l,k}^{z_{nl}^{w_{nk}}}

に比例します。w_{nk} は観測値です。

尤度関数を \beta 乗して  \phip_{l,k} にそれぞれハイパーパラメータ \alpha\gamma のディリクレ事前分布を置くと、ギブスサンプリングのアルゴリズムは以下のように導かれます。(まちがってたらすみません。)

 z_{nl} \sim \mathrm{Categorical}(\eta_n)
 \eta_{ln} = \{\phi_l p_{l,k}^{w_{nk}} / [\sum_{l=1}^{L} \phi_l p_{l,k}^{w_{nk}}] \} ^\beta

p_{l} \sim \mathrm{Dirichlet}(\sum_{n=1}^{N}\beta z_{n,l} w_{nk}+\alpha)

\phi \sim \mathrm{Dirichlet}(\sum_{n=1}^{N}\beta z_{n,l} +\gamma)

Rcpp で実装したコードは
multmix_Gibbs.cpp · GitHub
にあります。

実際にクラスタ数4の混合多項分布から乱数を生成して WBIC を計算してみると 4 で最小になりました。

f:id:abrahamcow:20190113171937p:plain

library(Rcpp)
sourceCpp("~/Documents/multmix_Gibbs.cpp")

set.seed(1)
prob1 <-  gtools::rdirichlet(4,rep(1,4))
clabel <- sample.int(4,100,replace = TRUE)
x <-t(sapply(clabel,function(cl)rmultinom(1,100,p=prob1[cl,])))

smplist <-lapply(2:9,function(l)multmix_Gibbs(x,l,rep(1/l,l),pini = gtools::rdirichlet(l,rep(1,4)),alpha = 1,iter = 2000,beta=1/log(100)))
lml <- sapply(smplist,function(x)-mean(x$loglik[-c(1:999)]))
plot(2:9,lml,xlab="L",ylab = "WBIC")
abline(h=min(lml),lty=2)
print((2:9)[which.min(lml)])