jagsで潜在混合モデル
ベイス統計で実戦モデリング(北大路書房)の第6章 潜在混合モデルの練習で、次のコードを作成してみた。
8個のデータが4個づつ違う正規分布に従うとして、データがどちらの正規分布に従うか、またその時の正規分布の平均を推測する。ただし、標準偏差はどちらも同じであるとする。
モデルは次のように作成した。
model1.txt
model {
for (i in 1: 8) {
z[i] ~ dbern(0.5)
}
lambda ~ dgamma(0.001, 0.001)
psi ~ dbeta(1,1)
m1 <- psi*0.1 + 0.48
phi ~ dbeta(1,1)
m2 <- phi*0.1 + 0.48
for (i in 1:8) {
mx[i] <- equals(z[i],0)*m1 + equals(z[i],1)*m2
k[i] ~ dnorm(mx[i], lambda)
}
}
zには1あるいは0が入り、どちらのグループかをこれで区別する。
以下がRのコードである。
rm(list=ls())
setwd("directory指定")
library(R2jags)
k <- c(30,26,11,45,28,10,9,8)
data <- list("k")
x <- c(rep(0,5000),rep(1,5000))
myinits <- list(list(z=sample(x, 8), psi <- 0.5, phi <- 0.5, lambda <- 0.1))
parameters <- c("z", "m1", "m2")
samples <- jags(data, inits=myinits, parameters,model.file="model1.txt", n.chains=1, n.burnin=1,n.iter=10000, n.thin=1, DIC=T)