Using an explicit likelihood...
model {
for (i in 1:N) {
Zero[i] <- 0
Zero[i] ~ dpois(phi[i])
phi[i] <- -log((p[1]*exp(-0.5*tau*pow(y[i] - lambda[1], 2))
+ p[2]*exp(-0.5*tau*pow(y[i] - lambda[2], 2)))
* sqrt(tau/(2*3.14159))) + C
}
C <- 100
p[2] <- 1 - p[1]
p[1] ~ dunif(0, 1)
theta ~ dunif(0, 1000)
lam ~ dunif(0, 1000)
lambda[2] <- lam + theta/2
lambda[1] <- lam - theta/2
sigma ~ dunif(0, 100)
tau <- 1/pow(sigma, 2)
dummy <- T[1] # so we can use the same data as Example 11.6.1 without WinBUGS complaining about T being undefined in the model
}
Note that this doesn't converge with standard Robert parameterisation - lambda[1] switches between two group means, and theta goes to infinity when lambda is in second group. Centred parameterisation given here works better, though mixing is slow.
Note also we cannot use p[1,2] ~ ddirch(). This would give the error "unable to choose update method" in WinBUGS, because the Dirichlet in WinBUGS can only be used as a conjugate prior for the categorical or multinomial, or for forward sampling.
Data:
list(y = c(529.0, 530.0, 532.0, 533.1, 533.4, 533.6, 533.7, 534.1, 534.8, 535.3,
535.4, 535.9, 536.1, 536.3, 536.4, 536.6, 537.0, 537.4, 537.5, 538.3,
538.5, 538.6, 539.4, 539.6, 540.4, 540.8, 542.0, 542.8, 543.0, 543.5,
543.8, 543.9, 545.3, 546.2, 548.8, 548.7, 548.9, 549.0, 549.4, 549.9,
550.6, 551.2, 551.4, 551.5, 551.6, 552.8, 552.9,553.2), N = 48,
T = c(1, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, 2))
Inits:
list(lam = 545, theta = 12, sigma = 20, p = c(0.5, NA))
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
Dbar Dhat pD DIC
Zero 9914.660 9911.010 3.656 9918.320
total 9914.660 9911.010 3.656 9918.320
Simple Normal model...
model {
for (i in 1:N) {
y[i] ~ dnorm(lam, tau)
}
lam ~ dunif(0, 1000)
sigma ~ dunif(0, 100)
tau <- 1/pow(sigma, 2)
dummy <- T[1]
}
Data:
list(y = c(529.0, 530.0, 532.0, 533.1, 533.4, 533.6, 533.7, 534.1, 534.8, 535.3,
535.4, 535.9, 536.1, 536.3, 536.4, 536.6, 537.0, 537.4, 537.5, 538.3,
538.5, 538.6, 539.4, 539.6, 540.4, 540.8, 542.0, 542.8, 543.0, 543.5,
543.8, 543.9, 545.3, 546.2, 548.8, 548.7, 548.9, 549.0, 549.4, 549.9,
550.6, 551.2, 551.4, 551.5, 551.6, 552.8, 552.9,553.2), N = 48,
T = c(1, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, 2))
Inits:
list(lam = 545, sigma = 20)
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
Dbar Dhat pD DIC
y 324.417 322.479 1.938 326.354
total 324.417 322.479 1.938 326.354
To compare to the mixture model (DIC=9918) add a constant 9600 = 2 * C * n = 2 * 100 * 48.