#  
# Data Summaries
# Quantities used in the likelihood
#
bary<-mean(bodytemp$temp)
n<-length(bodytemp$temp)
ss <- (n-1)*var(bodytemp$temp)
#
#
# ----------------
# sensitivity on c (prior mean =98.6
# ----------------

# prior parameters
prior<-list(  mu=98.6, c=100, a=0.001, b=0.001  )

prior$c <- 10^seq(-2, 5)

w                   <- n*prior$c/(n*prior$c+1) 
marginal.post.mu$mu <- w * bary + (1-w) * prior$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.tau$a <- 0.5*n+prior$a 
marginal.post.tau$b <- 0.5*ss + prior$b 

marginal.post.mu$sd <-  sqrt((w/n) * (ss + 2*prior$b)/(n+2*prior$a-2))


postq1 <- marginal.post.mu$mu + qt(0.025,marginal.post.mu$a) /sqrt(marginal.post.mu$tau)
postq3 <- marginal.post.mu$mu + qt(0.975,marginal.post.mu$a) /sqrt(marginal.post.mu$tau)


postscript('chap1_plot_sens1.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='r', bty='l' , cex=1.5)
plot( log(prior$c, 10) , marginal.post.mu$mu , type='l', ylab='Posterior Mean', xlab='Log10(c)', ylim=c(98.0,98.5))

lines(log(prior$c, 10) , postq1, lty=2)
lines(log(prior$c, 10) , postq3, lty=2)

graphics.off()
#
#
# ----------------
# sensitivity on c (prior mean =00.0
# ----------------

# prior parameters
prior<-list(  mu=0.0, c=100, a=0.001, b=0.001  )

prior$c <- 10^seq(-2, 5)

w                   <- n*prior$c/(n*prior$c+1) 
marginal.post.mu$mu <- w * bary + (1-w) * prior$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.tau$a <- 0.5*n+prior$a 
marginal.post.tau$b <- 0.5*ss + prior$b 

marginal.post.mu$sd <-  sqrt((w/n) * (ss + 2*prior$b)/(n+2*prior$a-2))


postq1 <- marginal.post.mu$mu + qt(0.025,marginal.post.mu$a) /sqrt(marginal.post.mu$tau)
postq3 <- marginal.post.mu$mu + qt(0.975,marginal.post.mu$a) /sqrt(marginal.post.mu$tau)


postscript('chap1_plot_sens1b.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
plot( log(prior$c, 10) , marginal.post.mu$mu , type='l', ylab='Posterior Mean', xlab='Log10(c)' , ylim=c(50,100))
#
#lines(log(prior$c, 10) , postq1, lty=2)
#lines(log(prior$c, 10) , postq3, lty=2)
graphics.off()

# ----------------
# sensitivity on mu
# ----------------


prior$c <- 100

prior$mu <-seq(98,99, 0.01)

w                   <- n*prior$c/(n*prior$c+1) 
marginal.post.mu$mu <- w * bary + (1-w) * prior$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.tau$a <- 0.5*n+prior$a 
marginal.post.tau$b <- 0.5*ss + prior$b 

marginal.post.mu$sd <-  sqrt((w/n) * (ss + 2*prior$b)/(n+2*prior$a-2))


postscript('chap1_plot_sens2.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
plot( prior$mu , marginal.post.mu$mu , type='l', ylab='Posterior Mean', xlab='Log10(c)', ylim=c(98.2,98.3))
graphics.off()