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)
    }
}

 

psiとphiが正規分布の平均値、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)