##############################################################################################################
#
#  Figure 3.2 (page 101) - Normal data with known population and Cauchy prior for population mean
#
#  Bayes' theorem via the rejection method in the context of Example 2.1 and Figure 2.1
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################

# Sufficient statistics
# ---------------------
xbar  = 7.0
sig2n = 4.5

# Hyperparameters of the normal prior
# -----------------------------------
mu0   = 0.0
tau20 = 1.0

# Hyperparameters of the conjugate normal posterior
# -------------------------------------------------
tau21 = 1/(1/sig2n+1/tau20)
mu1   = tau21*(xbar/sig2n+mu0/tau20)


# Sampling from the prior
# -----------------------
set.seed(21386)
M  = 200
x1 = mu0+sqrt(tau20)*rt(M,1)

# Figure 3.2
# ----------
x          = seq(-4,14,length=5000)
prior      = dt((x-mu0)/sqrt(tau20),1)/sqrt(tau20)
likelihood = dnorm(x,xbar,sqrt(sig2n))
posterior  = prior*likelihood
w          = dnorm(x1,xbar,sqrt(sig2n))
w          = w/sum(w)
plot(x,posterior,xlab="",ylab="",main="",type="l",axes=F,ylim=c(0,0.004))
axis(1,at=c(-4,0,5,10,14))
axis(2)
lines(x,likelihood/max(likelihood)*0.002,lty=2)
lines(x,prior/max(prior)*0.0036,lty=3)
for (i in 1:M){
  segments(x1[i],0,x1[i],-0.0001)
  segments(x1[i],0,x1[i],w[i]*0.005)
}
abline(h=0)