model {
   # genotype frequencies
   y[1] <- x[1]*x[1]
   y[2] <- 2*x[1]*x[2]
   y[3] <- 2*x[1]*x[3]
   y[4] <- 2*x[1]*x[4]
   y[5] <- x[2]*x[2]
   y[6] <- 2*x[2]*x[3]
   y[7] <- 2*x[2]*x[4]
   y[8] <- x[3]*x[3]
   y[9] <- 2*x[3]*x[4]
   y[10] <- x[4]*x[4]

   # gametic disequilibrium
   D.ab <- x[1]*x[4] - x[2]*x[3]

   # phenotype frequencies
   p[1] <- y[1]
   p[2] <- y[2]
   p[3] <- y[3]
   p[4] <- y[5]
   # the double heterozygotes
   p[5] <- y[4] + y[6]
   p[6] <- y[7]
   p[7] <- y[8]
   p[8] <- y[9]
   p[9] <- y[10]

   # likelihood
   n[1:9] ~ dmulti(p[], N)

   # prior
   # allele frequencies
   for (i in 1:4) {
      phi[i] ~ dgamma(1, 0.01)
   }
   for (i in 1:4) {
      x[i] <- phi[i]/sum(phi[])
   }

   # sample size
   N <- sum(n[])
}

list(n = c(91, 147, 85, 32, 78, 75, 5, 17, 7))