model {
    for(k in 1:K){
      beta[k] <- alpha[2,k]/sum(alpha[2,1:K]) - alpha[1,k]/sum(alpha[1,1:K])
    }

    for(i in 1:I1){
      ## reference group
      for(n in 1:N1[i]){
        xG1[n,1,i] ~ dcat(pStates1[1, 1:numStates, sCat1[i,n]])
        for(h in 1:numFlank){
          xG1[n,h+1,i] ~ dcat(pStates2[h, 1:4, sCat1[i,n]])
        }
        sCat1[i,n] ~ dcat(p1[i,1:K])
      }

      p1[i,1:K] ~ ddirch(alpha[1,1:K])

    }

    for(i in 1:I2){
      ## the other group
      for(n in 1:N2[i]){
        xG2[n,1,i] ~ dcat(pStates1[1, 1:numStates, sCat2[i,n]])
        for(h in 1:numFlank){
          xG2[n,h+1,i] ~ dcat(pStates2[h, 1:4, sCat2[i,n]])
        }
        sCat2[i,n] ~ dcat(p2[i,1:K])
      }

      p2[i,1:K] ~ ddirch(alpha[2,1:K])
    }

    ## priors
    for(k in 1:K){
      alpha[1,k] ~ dgamma(0.001,0.001)T(0.05,)
      alpha[2,k] ~ dgamma(0.001,0.001)T(0.05,)
    }

    ## signatures
    for(k in 1:K) {
      ## substituion type
      pStates1[1,1:numStates,k] ~ ddirch(gammas1[1:numStates])

      ## flanking bases
      for(w in 1:numFlank){
        pStates2[w,1:4,k] ~ ddirch(gammas2[1:4])
      }
    }

    ## priors
    for(q in 1:numStates){
        gammas1[q] <- 1
    }

    for(h in 1:4){
        gammas2[h] <- 1
    }
}
