model {
   ## genotype frequencies
   ##
   for (i in 1:n.pops) {
     for (j in 1:n.loci) {
       x[i,j,1] <- p[i,j]*p[i,j] + f*p[i,j]*(1-p[i,j])
       x[i,j,2] <- 2*(1-f)*p[i,j]*(1 - p[i,j])
       x[i,j,3] <- (1-p[i,j])*(1-p[i,j]) + f*p[i,j]*(1-p[i,j])
     }
   }

   ## likelihood
   ##
   for (i in 1:n.pops) {
     for (j in 1:n.loci) {
       n[i,j,1:3] ~ dmulti(x[i,j,], N[i,j])
     }
   }

   ## priors
   ##
   ## allele frequencies within populations
   ##
   for (i in 1:n.pops) {
     for (j in 1:n.loci) {
       p[i,j] ~ dbeta(alpha[j], beta[j])
     }
   }
   ## inbreeding coefficient within populations
   ##
   f ~ dunif(0, 1)

   ## theta (Fst)
   ##
   theta ~ dunif(0,1)

   ## pi
   ##
   for (i in 1:n.loci) {
     pi[i] ~ dunif(0,1)
   }

   ## parameters of the beta distribution
   ## the weird constraints are to ensure that both of them
   ## lie in [1, 1.0e4]
   ##
   for (i in 1:n.loci) {
     alpha[i] <- max(1, min(((1-theta)/theta)*pi[i], 1.0e4))
     beta[i]  <- max(1, min(((1-theta)/theta)*(1-pi[i]), 1.0e4))
   }
}