model {

  itodist[1,1,1,1:4] ~ ddirch(alphaito1[])
  itodist[1,1,2,1:4] ~ ddirch(alphaito2[])
  itodist[1,2,1,1:4] ~ ddirch(alphaito1[])
  itodist[1,2,2,1:4] ~ ddirch(alphaito2[])
  itodist[2,1,1,1:4] ~ ddirch(alphaito1[])
  itodist[2,1,2,1:4] ~ ddirch(alphaito2[])
  itodist[2,2,1,1:4] ~ ddirch(alphaito2[])
  itodist[2,2,2,1:4] ~ ddirch(alphaito2[])

  uetzdist[1,1,1:2] ~ ddirch(alphabin1[])
  uetzdist[1,2,1:2] ~ ddirch(alphabin1[])
  uetzdist[2,1,1:2] ~ ddirch(alphabin1[])
  uetzdist[2,2,1:2] ~ ddirch(alphabin2[])

  litdist[1,1,1:2] ~ ddirch(alphabin1[])
  litdist[1,2,1:2] ~ ddirch(alphabin1[])
  litdist[2,1,1:2] ~ ddirch(alphabin1[])
  litdist[2,2,1:2] ~ ddirch(alphabin2[])

  ppiprior[1:2] ~ ddirch(alphappi[])
  itofpprior[1:2] ~ ddirch(alphafp[])

  seeitoprior[1:2] ~ ddirch(alphasee[])
  seeuetzprior[1:2] ~ ddirch(alphasee[])
  seelitprior[1:2] ~ ddirch(alphasee[])

  for( p in 1 : N) {
       itopfp[p] ~ dcat(itofpprior[])
  }

  for( i in 1 : M ) {
       ppi[i] ~ dcat(ppiprior[])

       # Explaining away variable for ito,uetz,lit (if 0)
       seeito[i] ~ dcat(seeitoprior[])
       seeuetz[i] ~ dcat(seeuetzprior[])
       seelit[i] ~ dcat(seelitprior[])

       itofp[i] <- step(itopfp[parent1[i]] + itopfp[parent2[i]] - 3) + 1
       ito[i] ~ dcat(itodist[seeito[i],ppi[i],itofp[i],])
       uetz[i] ~ dcat(uetzdist[seeuetz[i],ppi[i],])
       lit[i] ~ dcat(litdist[seelit[i],ppi[i],])
  }

}