model
   {
      for (i in 1 : N){
         Y[i, 1 : 2] ~ dmnorm(mu[], tau[ , ])
      }
      mu[1] <- 0
      mu[2] <- 0
      tau[1 : 2,1 : 2] ~ dwish(R[ , ], 2)
      R[1, 1] <- 0.001 
      R[1, 2] <- 0
      R[2, 1] <- 0; 
      R[2, 2] <- 0.001
      Sigma2[1 : 2,1 : 2] <- inverse(tau[ , ])
      rho <- Sigma2[1, 2] / sqrt(Sigma2[1, 1] * Sigma2[2, 2])
   }