################################################################################
#
# NORMAL PRIOR AND NORMAL LIKELIHOOD LEADING TO NORMAL POSTERIOR
#
################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
################################################################################
# Prior hyperparameters for Physicists A and B
theta0A = 900
tau20A  = 20^2
theta0B = 800
tau20B  = 80^2
tau0A   = sqrt(tau20A)
tau0B   = sqrt(tau20B)

# Likelihood
sig2 = 40^2
x    = 850
sig  = sqrt(sig2)

# Posterior hyperparameters for Physicists A and B
tau21A  = 1/(1/tau20A+1/sig2)
tau21B  = 1/(1/tau20B+1/sig2)
theta1A = tau21A*(theta0A/tau20A+x/sig2)
theta1B = tau21B*(theta0B/tau20B+x/sig2)
tau1A = sqrt(tau21A)
tau1B = sqrt(tau21B)


# Densities
L = min(theta0A-4*tau0A,theta0B-4*tau0B,theta1A-4*tau1A,theta1B-4*tau1B)
U = max(theta0A+4*tau0A,theta0B+4*tau0B,theta1A+4*tau1A,theta1B+4*tau1B)
theta = seq(L,U,length=1000)

pA    = dnorm(theta,theta0A,tau0A)
pB    = dnorm(theta,theta0B,tau0B)
p     = dnorm(theta,x,sig)
poA   = dnorm(theta,theta1A,tau1A)
poB   = dnorm(theta,theta1B,tau1B)
U     = max(pA,pB,p,poA,poB)

#pdf(file="boxtiao.pdf",width=8,height=6)
plot(theta,pA,xlab=expression(theta),ylab="",type="l",ylim=c(0,U))
lines(theta,pA,col=2,lwd=2)
lines(theta,pB,col=3,lwd=2)
lines(theta,p,col=1,lwd=2)
lines(theta,poA,col=2,lwd=2,lty=2)
lines(theta,poB,col=3,lwd=2,lty=2)
legend(550,0.0175,legend=c("Physicist A: prior","Physicist B: prior","Likelihood","Physicist A: posterior","Physicist B: posterior"),col=c(2,3,1,2,3),lty=c(1,1,1,2,2),lwd=rep(2,5),bty="n")
#dev.off()