model {
  # convert data
  # N'goye
  n[1,1,1] <- idh1.1[1]
  n[1,1,2] <- idh1.1[2]
  n[1,1,3] <- idh1.1[3]
  n[2,1,1] <- idh2.1[1]
  n[2,1,2] <- idh2.1[2]
  n[2,1,3] <- idh2.1[3]
  n[3,1,1] <- pgm.1[1]
  n[3,1,2] <- pgm.1[2]
  n[3,1,3] <- pgm.1[3]
  n[4,1,1] <- mdh.1[1]
  n[4,1,2] <- mdh.1[2]
  n[4,1,3] <- mdh.1[3]
  n.sample[1] <- 50
  # Ghana
  n[1,2,1] <- idh1.2[1]
  n[1,2,2] <- idh1.2[2]
  n[1,2,3] <- idh1.2[3]
  n[2,2,1] <- idh2.2[1]
  n[2,2,2] <- idh2.2[2]
  n[2,2,3] <- idh2.2[3]
  n[3,2,1] <- pgm.2[1]
  n[3,2,2] <- pgm.2[2]
  n[3,2,3] <- pgm.2[3]
  n[4,2,1] <- mdh.2[1]
  n[4,2,2] <- mdh.2[2]
  n[4,2,3] <- mdh.2[3]
  n.sample[2] <- 40

  # likelihood
  # i - locus index
  # k - population index
  for (i in 1:4) {
    for (k in 1:2) {
       n[i,k,1:3] ~ dmulti(gamma[i,k,], n.sample[k])
    }
  }

  # genotype frequencies
  for (i in 1:4) {
    for (k in 1:2) {
      gamma[i,k,1] <- p[i,k]*p[i,k] + f*p[i,k]*(1-p[i,k])
      gamma[i,k,2] <- 2*p[i,k]*(1-p[i,k])*(1-f)
      gamma[i,k,3] <- (1-p[i,k])*(1-p[i,k]) + f*p[i,k]*(1-p[i,k])
    }
  }

  F <- 1 - (1-f)*(1-theta)

  # priors
  # within population allele frequencies
  for (i in 1:4) {
    alpha[i] <- ((1-theta)/theta)*pi[i]
    beta[i]  <- ((1-theta)/theta)*(1-pi[i])
    for (k in 1:2) {
      p[i,k] ~ dbeta(alpha[i], beta[i])
    }
  }
  # mean allele frequencies
  for (i in 1:4) {
    pi[i] ~ dunif(0, 1)
  }
  # f
  f ~ dunif(0, 1)
  # theta
  theta ~ dunif(0, 1)
}

# data
# pgm[3] = 33
# mdh[3] = 33
list(idh1.1=c(1, 5,44), idh2.1=c(0,2,48), pgm.1=c(32,18, 0), mdh.1=c(7,18,25),
     idh1.2=c(0,16,24), idh2.2=c(0,2,38), pgm.2=c(19,16, 5), mdh.2=c(0,13,27))

# inits
list(f=0.5, theta=0.5, pi=c(0.5,0.5,0.5,0.5))