model;
{
    mu ~ dnorm( 0.0,1.0E-6)
    rpsi ~ dgamma(0.001,0.001)    
    rphi ~ dgamma(0.001,0.001)
    psi <- 1/rpsi
    phi <- 1/rphi
    for( i in 1 : n ) {
      theta[i] ~ dnorm(mu,rpsi)
    }
   for (i in 1:N){
      x[i] ~ dnorm(theta[d[i]],rphi)
   }
}

data;
list(n=4,N=24,
x=c(62,60,63,59,
    63,67,71,64,65,66,
    68,66,71,67,68,68,
    56,62,60,61,63,64,63,59),
d=c(1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,4,4))

inits;
list(mu=60,rphi=0.1,rpsi=0.1,theta=c(60,60,60,60))