In R, I am running an MCMC Bayesian inference for data from mixture of Gamma distributions. JAGS is used here. The model file gmd.bug is as follows
model {
for (i in 1:N) {
y[i] ~ dsum(p*one, (1-p)*two)
}
one ~ dgamma(alpha1, beta1)
two ~ dgamma(alpha2, beta2) alpha1 ~ dunif(0, 10)
beta1 ~ dunif(0, 10)
alpha2 ~ dunif(0, 10)
beta2 ~ dunif(0, 10)
p ~ dunif(0, 1)
}
Then, this is the inference phase
gmd.jags = jags.model("gmd.bug",
data = list(N = NROW(a), y=unlist(a)),
inits = inits, n.chains = 3, n.adapt = 1000)
Here is the error that puzzled me
Error in jags.model("gmd.bug", data = list(N = NROW(a), y = unlist(a)), :
Error in node y[1]
Node inconsistent with parents
Anyone know what needs to be fixed here?
Answer to OP's original question
When you write y[i] ~ dsum(p*dgamma(alpha1, beta1), (1-p)*dgamma(alpha2, beta2))
, dgamma(alpha1, beta1)
needs to be indexed by [i], as in
gamma1[i] ~ dgamma(alpha1, beta1)
gamma2[i] ~ dgamma(alpha2, beta2)
Answer to OP's second question (after edits)
This is the crux of your problem. But fixing it raises additional difficulties, because in order to ensure that y[i] is consistent with its parents at initialization, you need to make sure that at initialization it is strictly true that y[i] == p*gamma1[i]+(1-p)*gamma2[i]
. If you let JAGS handle the initialization automatically, it will initialize from the priors, without understanding the constraint on initial values imposed by dsum
, and you will get an error. One strategy is to initialize both gamma1
and gamma2
at y
.
The following code works for me (but of course you'll want to run many more iterations):
# Data simulation:
library(rjags)
N=200
alpha1 <- 3
beta1 <- 3
alpha2 <- 5
beta2 <- 1
p <- .7
y <- vector(mode="numeric", length=N)
for(i in 1:N){
y[i] <- p*rgamma(1,alpha1,beta1) + (1-p)*rgamma(1,alpha1,beta1)
}
# JAGS model
sink("mymodel.txt")
cat("model{
for (i in 1:N) {
gamma1[i] ~ dgamma(alpha1, beta1)
gamma2[i] ~ dgamma(alpha2, beta2)
pg1[i] <- p*gamma1[i]
pg2[i] <- (1-p)*gamma2[i]
y[i] ~ dsum(pg1[i], pg2[i])
}
alpha1 ~ dunif(0, 10)
beta1 ~ dunif(0, 10)
alpha2 ~ dunif(0, 10)
beta2 ~ dunif(0, 10)
p ~ dunif(0, 1)
}", fill=TRUE)
sink()
jags.data <- list(N=N, y=y)
inits <- function(){list(gamma1=y, gamma2=y)}
params <- c("alpha1", "beta1", "alpha2", "beta2", "p")
nc <- 5
n.adapt <-200
n.burn <- 200
n.iter <- 1000
thin <- 10
mymodel <- jags.model('mymodel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt)
update(mymodel, n.burn)
mymodel_samples <- coda.samples(mymodel,params,n.iter=n.iter, thin=thin)