title: Hierarchical Mixture Model in Rjags
body: <p>This is my first time using rjags, and I'm trying to fit some count data Y. I'm using a hierarchical mixture model as follows:</p>
  
  <pre><code>Y ~ p*Poisson(N*lambda1) + (1-p)*Poisson(N*lambda2)
  lambda1 ~ Gamma(a,b)
  lambda2 ~ Lognormal(c,d)
  a ~ Gamma(1,1)
  b ~ Gamma(1,1)
  c ~ Normal(0,1)
  d ~ Gamma(1,1)
  </code></pre>
  
  <p>Here, Y is the count data that I observe, and N is known. </p>
  
  <p>I wrote up a simple rjags model that I've been playing around with. However, when testing on simple simulated data, I'm getting really bad results. Here is my code to generate the simulated data and run the model:</p>
  
  <pre><code>a &lt;- 0.5
  b &lt;- 0.5
  c &lt;- -10
  d &lt;- 1
  
  lambda1 &lt;- rgamma(30,a,b)
  lambda2 &lt;- rlnorm(70,c,d)
  counts &lt;- rpois(100,1000*c(lambda1,lambda2))
  
  model_string &lt;- "model{
    # Likelihood 
    for (i in 1:n) {
      mu1[i] &lt;- N*lambda1[i]
      mu2[i] &lt;- N*lambda2[i]
      lambda1[i] ~ dgamma(a,b)
      lambda2[i] ~ dlnorm(c,d)
      m[i] ~ dcat(mprior[])
      mu[i] &lt;- equals(m[i],1)*mu1[i] + equals(m[i],2)*mu2[i]
      Y[i] ~ dpois(mu[i])
    }
    # Prior
    mprior[1] &lt;- 0.5
    mprior[2] &lt;- 0.5
    a ~ dgamma(1,1)
    b ~ dgamma(1,1)
    c ~ dnorm(0,1)
    d ~ dgamma(1,1)
  }"
  
  model &lt;- jags.model(textConnection(model_string), 
                      data = list(Y=counts,N=1000,n=100))
  update(model,10000)
  samp &lt;- coda.samples(model, 
                       variable.names=c("a","b","c","d","m"), 
                       n.iter=20000)
  
  print(colMeans(samp[[1]])[1:4])
  </code></pre>
  
  <p>After running this, the posterior estimates of a, b, c, d are quite poor, and the component assignments m also don't match up well with the true assignments. I also notice that plotting the chains don't look great, even when increasing the number of iterations. </p>
  
  <p>Any advice? I'm not sure if I'm fitting the mixture in the best way. I'm also definitely open to changing my priors on a, b, c, d if there are other distributions that will be easier to work with.</p>