postplot <- function(mu, tau, ybar, sigma, n){
  denom <- sigma^2/n + tau^2
  mu <- (sigma^2/n*mu + tau^2*ybar)/denom
  stdev <- sqrt(sigma^2/n*tau^2/denom)
  return(c(mu,stdev))
  }

x <- seq(-2,8,length.out=100)
plot(x, dnorm(x, mean=2, sd=1), ylim=c(0, 1.3),
  xlab=expression(theta), ylab="density", type= "l")

param1 <- postplot(2,1,6,1,1)
lines(x, dnorm(x, mean=param1[1], sd=param1[2]), lty=2)
param10 <- postplot(2,1,6,1,10)
lines(x, dnorm(x, mean=param10[1], sd=param10[2]), lty=3)
legend(-2, 1.3, legend=c("prior", "posterior with n=1",
  "posterior with n=10"), lty=1:3, ncol=1)

y1 <- rnorm(2000, mean=param1[1], sd=param1[2])
## param1: [1] 4.0000000 0.7071068
r1 <- hist(y1, freq=F, breaks=20, plot=F)
lines(r1, lty=2, freq=F, col="gray90")