混合ポアソン分布のための崩壊型ギブスサンプリング

PoissonMixtureCollapsedGibbs.R

#このソースコードの自由な複製・改変・再配布を認める。
set.seed(1234)
N <- 100
ind <- sample.int(2,N,replace = TRUE)
lambda <- c(2,10)
y <- rpois(N,lambda[ind])
a <- 1 #ガンマ事前分布のパラメータ
b <- 1 #ガンマ事前分布のパラメータ
alpha <-c(2,2) #ディリクレ事前分布のパラメータ
softmax <- function(x){
  maxx <- max(x)
  exp(x-maxx)/sum(exp(x-maxx))
}
maxit <- 1000
shist <- array(NA_integer_,dim = c(N,2,maxit))
lphist <- numeric(maxit-1)
s <- t(rmultinom(N,1,c(0.5,0.5))) #潜在変数の初期化
alphahat <- colSums(s)+alpha
ahat <- colSums(s*y)+a
bhat <- colSums(s)+b
shist[,,1]<-s
for(i in 2:maxit){
  for(n in 1:N){
    ahat <- ahat - y[n]*shist[n,,i-1]
    alphahat <- alphahat - shist[n,,i-1]
    bhat <- bhat - shist[n,,i-1]
    unnorm <-dnbinom(y[n],ahat,1-1/(1+bhat),log = TRUE)+log(alphahat)
    lphist[i-1] <- lphist[i-1] + log(sum(exp(unnorm)))
    snew <- drop(rmultinom(1,1,softmax(unnorm)))
    shist[n,,i] <- snew
    alphahat <- alphahat + snew 
    ahat <- ahat + snew*y[n]
    bhat <- bhat + snew
  }
}
plot(lphist,type="l")
shat <-apply(shist[,,-c(1:100)], 1:2, mean)
cluster <-factor(apply(shat, 1, which.max))
print(mean(ind==cluster))

library(ggplot2)
df <- data.frame(y,cluster)
p1 <- ggplot(df,aes(x=y,fill=cluster,colour=cluster))+
  geom_histogram(position = "identity",alpha=0.5,bins=15)
print(p1)