# comments
x = c(seq(39.8 - 4*2, 39.8, length=100), seq(39.8, 39.8 + 4*2, length=100))
#generate a sequence to plot the density
#note that the mode is explicitly included twice in the sequence

# plot the prior distribution
plot(x, dnorm(x, mean=38, sd=sqrt(9)), type="l",  main="Prior Density", ylab="density")
# dnorm gives the normal density
# type  = l connects points by lines

# plot the standardized likelihood
plot(x, dnorm(x, mean=39.8, sd=sqrt(4/10000)), type="l", ylab="",main="Likelihood")
lines(x, dnorm(x, mean=38, sd=sqrt(9)), lty=2, col=2)
# add the prior density
lines(x, dnorm(x, mean=39.8, sd=sqrt(1/2500.111)), col=3)
# overlay the posterior density
x = seq(39.7, 39.9, length=100)

#Plot the posterior
plot(x, dnorm(x,mean=39.8, sd=sqrt(1/2500.111)), type="l", ylab="density",main="Posterior Density", xlab=expression(mu))
     
# Calculate an HPD Region
hl = qnorm(.025, 39.8, sqrt(1/2500.111))
hu = qnorm(.975, 39.8, sqrt(1/2500.111))
#> hu 
#[1] 39.8392
#> hl
#[1] 39.7608
 
hd = dnorm(hu, 39.8, sqrt(1/2500.111))
hd = dnorm(hl, 39.8, sqrt(1/2500.111))
#add a horzontal line at the height of highest densities
abline(h= hd)
# add verticle line segments at the points of highest densities
segments(c(hl, hu), c(0, 0), c(hl, hu), c(hd, hd))