#prior#
sd0 <- 3
phi0 <- sd0^2
theta0 <- 38
z=seq(28,47,by=0.05)
y0=dnorm(z,theta0,sd0)
plot(z,y0,type="l",ylim=c(0,0.9))

#likelihood#
sdobs <- 2
phi <- sdobs^2
nobs <- 10000     #NB: misprint here in Lee's code (chests.txt)# 
xbar <- 39.8

#posterior#
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")
y1=dnorm(z,theta1,sd1)
lines(z,y1)

#predictive distribution#
sd2=sqrt(phi+phi1)
cat("Predictive mean",theta1,"and s.d.",sd2,".\n")
y2=dnorm(z,theta1,sd2)
lines(z,y2)