Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
混合ポアソン分布のための崩壊型ギブスサンプリング
#このソースコードの自由な複製・改変・再配布を認める。
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment