model;
{
   for( k in 1 : K ) {
      for( m in 1 : 4 ) {
         logit(tau[k , m]) <- gamma[k , m]
         gamma[k , m] ~ dnorm( 0.0,0.45)
      } 
   }

   alpha[1] ~ dnorm(0,0.45)
   alpha[2] ~ dnorm(0,0.45)
   alpha[3] ~ dnorm(0,0.45)
   pi[1] <- exp(alpha[1]) / (1 + exp(alpha[1]) + exp(alpha[2]) + exp(alpha[3]))
   pi[2] <- exp(alpha[2]) / (1 + exp(alpha[1]) + exp(alpha[2]) + exp(alpha[3]))
   pi[3] <- exp(alpha[3]) / (1 + exp(alpha[1]) + exp(alpha[2]) + exp(alpha[3]))
   pi[4] <- 1 - pi[1] - pi[2] - pi[3] 

for(i in 1:N) {
	eta[i] ~ dcat(pi[1:4])
	for(k in 1:K) {
	    logit(p[i , k]) <- gamma[k , eta[i]]
         y[k,i] ~ dbin(p[i , k],1)
      }

   }
}