CODEWITHSUNDEEP
rmcmcjags
lilkaskitc
lilkaskitc

Reputation: 15

JAGS error for MCMC Bayesian inference

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?

Upvotes: 0

Views: 1244

Answers (1)

Jacob Socolar
Jacob Socolar

Reputation: 1202

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)

Upvotes: 1

Related Questions