sd0 <- 3
phi0 <- sd0^2
theta0 <- 38
sdobs <- 2
phi <- sdobs^2
nobs <- 1000
xbar <- 39.8
phi1 <- (phi0^(-1) + (phi/nobs)^(-1))^(-1)
sd1 <- sqrt(phi1)
theta1 <- phi1*(theta0/phi0 + xbar/(phi/nobs))
cat("Posterior mean",theta1,"and s.d.",sd1,".\n")