model {
   # Priors
   j ~ dcat(p[])
   psi ~ dbeta(1,1)
   rho ~ dbeta(100,100)
   for (i in 1:3){
    p[i] <- 0.33
    theta[i] <- ifelse(i==j, psi, rho) 
    k[i] ~ dbin(theta[i],n)
   }
   
}