model {

  mu ~ dnorm( 0.0, 1.0E-6 )
  tau ~ dgamma( 0.001, 0.001 )
 
  for ( i in 1:n ) {

    y[ i ] ~ dnorm( mu, tau )
 
  }

  sigma <- 1.0 / sqrt( tau )
  y.new ~ dnorm( mu, tau )
 
}