model {
   # allele frequencies
   p.a <- a[1] + a[2]/2.0
   p.b <- b[1] + b[2]/2.0

   # one-locus genotype frequencies
   a[1] <- p[1] + p[2] + p[3]
   a[2] <- p[4] + p[5] + p[6]
   a[3] <- p[7] + p[8] + p[9]
   b[1] <- p[1] + p[4] + p[7]
   b[2] <- p[2] + p[5] + p[8]
   b[3] <- p[3] + p[6] + p[9]

   # one-locus disequilibria
   D.a <- a[1] - p.a*p.a
   D.b <- b[1] - p.b*p.b
   f.a <- D.a/(p.a*(1-p.a))
   f.b <- D.b/(p.b*(1-p.b))

   # composite digenic disequilibrium
   Delta.ab <- 2*p[1] + p[2] + p[4] + p[5]/2 - 2*p.a*p.b

   # trigenic disequilibrium
   p.aab <- p[1] + p[2]/2
   p.abb <- p[1] + p[4]/2
   D.aab <- p.aab - p.a*Delta.ab - p.b*D.a - p.a*p.a*p.b
   D.abb <- p.abb - p.b*Delta.ab - p.a*D.b - p.a*p.b*p.b

   # quadrigenic disequilbrium
   Delta.aabb <- p[1] - 2*p.a*D.abb - 2*p.b*D.aab - 2*p.a*p.b*Delta.ab - Delta.ab*Delta.ab - p.a*p.a*D.b - p.b*p.b*D.a - D.a*D.b - p.a*p.a*p.b*p.b

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

   # prior
   # phenotype frequencies
   for (i in 1:9) {
      phi[i] ~ dgamma(1, 1)
   }
   for (i in 1:9) {
      p[i] <- phi[i]/sum(phi[])
   }

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

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