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, beta)
     }
   }
   ## inbreeding coefficient within populations
   ##
   f <- 0.0

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

   ## pi
   ##
   pi ~ dunif(0,1)

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