#WinBUGS model:
model{
   for(i in 1:N){Y[i,1:M] ~ dmnorm(mu[],InvCov[,])}                                  
   for(j in 1:M){mu[j] <- b[1]+b[2]*age[j]}
   for(j in 1:2){b[j] ~ dnorm(0.0, 0.001)}

    #Covariance models:
     for(j1 in 1:M){for(j2 in 1:M){
       Cor[j1,j2,1]<-equals(j1,j2)                #Ind
       Cor[j1,j2,2]<-rho1+equals(j1,j2)*(1-rho1)  #CS
       Cor[j1,j2,3]<-pow(rho2,abs(j1-j2))         #AR
       Cov[j1,j2]<-pow(sig,2)*Cor[j1,j2,type]
     }}
     InvCov[1:M,1:M]<-inverse(Cov[,])

    sig~dunif(0,10)
    rho1~dbeta(1,1)
    rho2~dbeta(1,1)
    type<-3                   #type=3 -> AR
    # type~dcat(p[])          #Stochastic search
    # for(j in 1:3){p[j]<-1/3}

    #Wishart prior:
#   InvCov[1:M,1:M]  ~ dwish(R[,], 4)	
#   for(j1 in 1:M){for(j2 in 1:M){R[j1,j2]<-equals(j1,j2)}}
#    Cov[1:M,1:M] <- inverse(InvCov[,])

   #Keep track of the correlation:
#   for(j1 in 1:M){for(j2 in 1:M){
#          Cor[j1,j2]<-Cov[j1,j2]/sqrt(Cov[j1,j1]*Cov[j2,j2])
#   }}

}

#Initial value
list(type=1)

#Jaw data
list(M = 4, N = 20, age = c(8.0, 8.5, 9.0, 9.5),
     Y = structure(.Data = c(47.8, 48.8, 49.0, 49.7,
	      	             46.4, 47.3,  47.7, 48.4,
			     46.3, 46.8, 47.8, 48.5,
			     45.1, 45.3, 46.1, 47.2,
			     47.6, 48.5, 48.9, 49.3,
			     52.5, 53.2, 53.3, 53.7,
			     51.2, 53.0, 54.3, 54.5,
			     49.8, 50.0, 50.3, 52.7,
			     48.1, 50.8, 52.3, 54.4,
			     45.0, 47.0, 47.3, 48.3,
			     51.2, 51.4, 51.6, 51.9,
			     48.5, 49.2, 53.0, 55.5,
			     52.1, 52.8, 53.7, 55.0,
			     48.2, 48.9, 49.3, 49.8,
			     49.6, 50.4, 51.2, 51.8,
			     50.7, 51.7, 52.7, 53.3,
			     47.2, 47.7, 48.4, 49.5,
			     53.3, 54.6, 55.1, 55.3,
			     46.2, 47.5, 48.1, 48.4,
			     46.3, 47.6, 51.3, 51.8),.Dim = c(20, 4))