model; { for( i in 1 : m ) { x[i] ~ dnorm(lambda,psi) } lambda ~ dnorm( 0.0,1.0E-6) mu ~ dnorm( 0.0,1.0E-6) psi ~ dgamma(0.01,0.001) omega ~ dgamma(0.01,0.001) kappa <- omega/psi for( j in 1 : n ) { y[j] ~ dnorm(mu,omega) } } data; list(m=12,n=8,x=c(35,26,24,12,27,17,-14,10,1,24,10,28), y=c(343,90,16,382,69,140 ,49 ,89)) inits; list(lambda=0,mu=0,psi=0.1,omega=0.1) (2.310 subtracted from A data and 2.298 from C data, and both sets multiplied by 100000 so as not to lose precision).