model{ 
      y[1:s] ~ dmulti(th[1 : s] , n)
      sum.g <- sum(g[])
    # smoothed frequencies
       for (i in 1 : s) { 
         Sm[i] <- n * th[i]
          g[i] <- exp(gam[i]) 
         th[i] <- g[i] / sum.g
      }
    # prior on elements of AR Precision Matrix 
       rho ~ dunif(0, 1)
       tau ~ dunif(0.5, 10)
    # MVN for logit parameters
       gam[1 : s] ~ dmnorm(mu[], T[ , ])
       for (j in 1:s) { 
         mu[j] <- -log(s)
      }
   # Define Precision Matrix
      for (j in 2 : s - 1) {
         T[j, j] <- tau * (1 + pow(rho, 2))
      }
      T[1, 1] <- tau 
      T[s, s] <- tau
      for (j in 1 : s -1 ) { 
         T[j, j + 1] <- -tau * rho
         T[j + 1, j] <- T[j, j + 1]
      }
      for (i in 1 : s - 1) {
         for (j in 2 + i : s) {
            T[i, j] <- 0; T[j, i] <- 0 
         }
      }
    # Or Could do in terms of covariance, which is simpler to write but slower
    #      for (i in 1 : s) {
    #         for (j in 1 : s) {
    #            cov[i, j] <- pow(rho, abs(i - j)) / tau
    #         }
    #      }
    #      T[1 : s, 1 : s] <- inverse(cov[ , ])
   }