model{ alpha <- mu^2/stdev^2 beta <- stdev^2/mu vartheta <- rphi/n0 theta ~ dnorm(theta0,vartheta) rphi ~ dgamma(alpha,beta) phi <- 1/rphi for (i in 1:n){ x[i] ~ dnorm(theta,rphi) } } data; list(n=12,n0=15,theta0=110,mu=300,stdev=160, x=c(141,102,73,171,137,91,81,157,146,69,121,134))