
library(rjags)

model = "tmp_model.bug"  
write("
model {
   for (i in 1:length(x)) {
      x[i] ~ dnorm(mu.s, tau)
   }
   mu.s <- mu + z;
   mu ~ dnorm(0.0, 1.0E-6);
   z ~ dnorm(0, 1/sigma.z^2)
   tau ~ dgamma(1.0, 1.0E-4)
   sigma <- 1/sqrt(tau)
   x.f ~ dnorm(mu.s, tau)
} 
", model)

mu.true = 3
sigma.true = 2
n = 20
x <- rnorm(n, mu.true, sigma.true)

sigma.z = 0.5

data <- list(x=x, sigma.z=sigma.z)

jm <- jags.model(model, data)
update(jm, 100)
chain <- coda.samples(jm, c('mu', 'sigma', 'x.f'), n.iter=10000)


plot(chain)
print(summary(chain))

cat(sprintf("mean(x) = %.2f,  sd(x) = %.2f, sd(x)/sqrt(n) =  %.2f\n",
            mean(x), sd(x), sd(x)/sqrt(n)))
