#  
# Data Summaries
# Quantities used in the likelihood
#
bary<-mean(bodytemp$temp)
n<-length(bodytemp$temp)
ss <- (n-1)*var(bodytemp$temp)
#
#
# prior parameters
prior<-list(  mu=98.6, c=0.01, a=1, b=1  )

postscript('chap1_ex4_inf1.ps', horizontal=FALSE, width = 12.0, height = 8.0,)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5, mfrow=c(2,2))

cvalues<-c(0.1,0.01,0.001,0.0001)

for (i in cvalues) {

prior$c<-i

#
#
# full posterior parameters
#
full.post<-list()

w <- n*prior$c/(n*prior$c+1)

full.post$mu <- w * bary + (1-w) * prior$mu 
full.post$c  <- w/n
full.post$a  <- prior$a + n/2
full.post$b  <- prior$b + ss/2
# ---------------------------------------
#
# Marginal posterior parameters for mu
# ---------------------------------------
marginal.post.mu<-list()
marginal.post.mu$mu <- full.post$mu
marginal.post.mu$tau   <- (n/w) * (n+2*prior$a)/(2*prior$b + ss)
marginal.post.mu$a <- n+2*prior$a 

marginal.post.mu$sd <-  sqrt((w/n) * (ss + 2*prior$b)/(n+2*prior$a-2))
#
# Plot Marginal posterior density
#
xrange<-5

prior$tau <-  (prior$a / prior$b )/prior$c


x<-seq(98.0,98.8, 0.001)

z <- (x - marginal.post.mu$mu) * sqrt(marginal.post.mu$tau) 
y <- dt( z,  marginal.post.mu$a ) * sqrt(marginal.post.mu$tau) 
# y <- ( 1+ (marginal.post.mu$tau/marginal.post.mu$a) * (x-marginal.post.mu$mu)^2 )^(-(marginal.post.mu$a+1)/2 )



zprior <- (x - prior$mu) * sqrt(prior$tau) 
yprior <- dt( zprior,  2*prior$a ) * sqrt(prior$tau) 
#yprior <- ( 1+ (prior$tau/prior$a) * (x-prior$mu)^2 )^(-(prior$a+1)/2 )

#y3<-dnorm( x, prior$mu, 1/sqrt(prior$tau) )
#y3<-dnorm( zprior, 0, 1)* sqrt(prior$tau) 

ydata <- dnorm( x, bary, sqrt( var(bodytemp$temp)/n ) )

plot(x,y,type='l', xlab='Mean of Temperature', ylab='Density Values', xlim=range(x), ylim=range( y, yprior ) )
title( paste('c=', prior$c))
lines(x,yprior, col=2, lty=2)

lines(x,ydata, col=3, lty=4)
}

graphics.off()