model {

   p[1:numPrey] ~ ddirch(alpha[]);  # these are weights for means
   for(i in 1:numPrey) {
      p2[i] <- p[i]*p[i]; # these are weights for variances
   }

   for(i in 1:numIsotopes) {
      # uniform prior on sds
      #sigma[i] ~ dunif(0,50);
      #isigma2[i] <- 1/(sigma[i]*sigma[i]);

      # inv gamma prior on variances
      #isigma2[i] ~ dgamma(0.01,0.01)
      #sigma[i] <- sqrt(1/isigma2[i]);

      # half-cauchy on sds
      tau.sigma[i] ~ dgamma(0.5,0.5);
      xi.sigma[i] ~ dnorm(0,1)
      sigma[i] <- abs(xi.sigma[i])/sqrt(tau.sigma[i]);
      isigma2[i] <- 1/(sigma[i]*sigma[i]);
      resid.var[i] <- sigma[i]*sigma[i];
   }

   # for each isotope, calculate the predicted mixtures
   for(i in 1:numIsotopes) {
      mix.mu[i] <- inprod(u[,i],p[]);
      mix.var[i] <- inprod(vars[,i],p2[]);
      mix.totalVar[i] <- mix.var[i] + resid.var[i];
      mix.prcsn[i] <- 1/mix.totalVar[i];
   }
   # This section does the likelihood / posterior, looping over N data points
   for(i in 1:N) {
      for(j in 1:numIsotopes) {
         X[i,j] ~ dnorm(mix.mu[j], mix.prcsn[j]);
      }
   }

}