model {

   # This section is the priors for the dirichlet...this is exactly the same as MixSIR
   # region-specific means
   p[1:num.prey] ~ ddirch(alpha[]);  
   for(i in 1:num.prey) {
      p2[i] <- p[i]*p[i]; # these are weights for variances
   }

   # v.n = v.0 + n, sig2|y ~ Inv-X2(v.n, s2)
   # to simulate, draw X from X2(v.n), theta=v.n*s2/X
   # the noninformative prior has v.0 = 0
   for(prey in 1:num.prey) {
   for(iso in 1:num.iso) {
      # this is the variance parameter
      # described Gelman p480 ("Inverse Chi-Square")
      temp.X[prey,iso] ~ dchisqr(sigma.v0.n[prey,iso]);
      bayes.sig2[prey,iso] <- sigma.v0.n[prey,iso]*sigma2.n[prey,iso]/temp.X[prey,iso];
      # calculate precision as reciprocal of variance
      bayes.iSig2[prey,iso] <- 1/bayes.sig2[prey,iso];
            
      # this is the scale parameter...below 3.7, Gelman p72
      k.n[prey,iso] <- u.scale[prey,iso] + samp.size[prey,iso];
      # this is just the equation for u.n, u below is the observed sample mean
      u.n[prey,iso] <- (u.scale[prey,iso]/k.n[prey,iso])*u.loc[prey,iso] + (samp.size[prey,iso]/k.n[prey,iso])*u[prey,iso];
      # precision param for mean
      temp.isig2[prey,iso] <- k.n[prey,iso]*bayes.iSig2[prey,iso];
      # this is equation 3.8, Gelman p72
      bayes.mu[prey,iso] ~ dnorm(u.n[prey,iso], temp.isig2[prey,iso]);
   }
   }

   # include variation in the mean and variance
   for(iso in 1:num.iso) {
      mix.mu[iso] <- inprod(bayes.mu[1:num.prey,iso],p[1:num.prey]);
      mix.var[iso] <- inprod(bayes.sig2[1:num.prey,iso],p2[1:num.prey]);
      mix.prcsn[iso] <- 1/(mix.var[iso]);
   }

   # This section does the likelihood / posterior, N data points
   for(i in 1:N) {
      for(iso in 1:num.iso) {
         X[i,iso] ~ dnorm(mix.mu[iso], mix.prcsn[iso]);
      }
   }

}