model {
   # genotype frequencies
   y[1] <- x[1]*x[1] + f*x[1]*(1-x[1])
   y[2] <- 2*x[1]*x[2]*(1-f)
   y[3] <- 2*x[1]*x[3]*(1-f)
   y[4] <- 2*x[1]*x[4]*(1-f)
   y[5] <- x[2]*x[2] + f*x[2]*(1-x[2])
   y[6] <- 2*x[2]*x[3]*(1-f)
   y[7] <- 2*x[2]*x[4]*(1-f)
   y[8] <- x[3]*x[3] + f*x[3]*(1-x[3])
   y[9] <- 2*x[3]*x[4]*(1-f)
   y[10] <- x[4]*x[4] + f*x[4]*(1-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, 1)
   }
   for (i in 1:4) {
      x[i] <- phi[i]/sum(phi[])
   }
   # f
   w ~ dunif(0,1)
   for (i in 1:4) {
     x.f[i] <- max(-x[i]/(1-x[i]), -(1-x[i])/x[i])
   }
   min.f.1 <- max(x.f[1], x.f[2])
   min.f.2 <- max(x.f[3], x.f[4])
   f.min <- max(min.f.1, min.f.2)
   f <- w*(1 - f.min) + f.min
   

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

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