model;
{
   mu[1:2] ~ dmnorm(mu0[],T[,])
   T[1:2,1:2] ~ dwish(R[,],2)
   Sigma[1:2,1:2] <- inverse(T[,])
  rho <- Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])
   for (i in 1:n){
      x[i,1:2] ~ dmnorm(mu[],T[,])
   }
}

data;
list(n=9,x=structure(.Data=c(22.5,17.0,20.1,14.9,23.3,16.0,22.9,17.4,23.1,
                             17.4,22.0,16.5,22.3,17.2,23.6,17.2,24.7,18.0),
                     .Dim=c(9,2)),
     mu0=c(20,15),
     R = structure(.Data = c(1,0,0,1),.Dim = c(2,2)))