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];
}
# these are for the mixture of variances
mixP[1:2] ~ ddirch(alpha.mix[]);
G ~ dcat(mixP[])
G.01 <- G-1
# 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] <- (1-G.01)*mix.var[i] + (G.01)*(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]);
}
}
}