#  
# Data Summaries
# Quantities used in the likelihood
#
source('chap01/ex5/chap01_ex5_bodytemp_data.r')
bodytemp<-ex1_5.bodytemp
bary<-mean(bodytemp$temp)
n<-length(bodytemp$temp)
ss <- (n-1)*var(bodytemp$temp)
#
#
# prior parameters
prior<-list(  mu=98.6, c=100, a=0.001, b=0.001  )
#
#
# 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

x <- seq( marginal.post.mu$mu - xrange*marginal.post.mu$sd, marginal.post.mu$mu + xrange*marginal.post.mu$sd, 0.1*marginal.post.mu$sd )
y <- ( 1+ (marginal.post.mu$tau/marginal.post.mu$a) * (x-marginal.post.mu$mu)^2 )^(-(marginal.post.mu$a+1)/2 )


par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
plot(x,y,type='l', xlab='Mean of Temperature', ylab='Unormalized posterior density')

postscript('chap1_plot1.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
plot(x,y,type='l', xlab='Mean of Temperature', ylab='Unormalized posterior density')

graphics.off()

# ---------------------------------------
#
# Marginal posterior parameters for tau
# ---------------------------------------

marginal.post.tau<-list()
marginal.post.tau$a <- 0.5*n+prior$a 
marginal.post.tau$b <- 0.5*ss + prior$b 

marginal.post.tau$mean <- marginal.post.tau$a/marginal.post.tau$b
marginal.post.tau$sd   <- sqrt(marginal.post.tau$a)/marginal.post.tau$b
#
#
# Plot of Marginal posterior density for tau
#
#
x<-seq( max(c(0, marginal.post.tau$mean - xrange*marginal.post.tau$sd)), marginal.post.tau$mean + xrange*marginal.post.tau$sd, 0.1*marginal.post.tau$sd)
y<- dgamma( x, marginal.post.tau$a, marginal.post.tau$b )

plot(x,y,type='l', xlab='Precision of Temperature', ylab='Posterior density')

postscript('chap1_plot2.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
plot(x,y,type='l', xlab='Precision of Temperature', ylab='Posterior density')
graphics.off()

# ---------------------------------------
#
# Marginal posterior parameters for s2
# ---------------------------------------

marginal.post.s2<-marginal.post.tau

marginal.post.s2$mean <- (ss/2 + prior$b)/( n/2 + prior$a -1) 
marginal.post.s2$sd   <- sqrt((ss/2 + prior$b)^2 /((n/2 + prior$a-1)^2 *(n/2+prior$a-2) ))
#
#
# Plot of Marginal posterior density for s2
#
#
x<-seq( max(c(0, marginal.post.s2$mean - xrange*marginal.post.s2$sd)), marginal.post.s2$mean + xrange*marginal.post.s2$sd, 0.1*marginal.post.s2$sd)

logconst<- marginal.post.s2$a * log(marginal.post.s2$b) - lgamma(marginal.post.s2$a)
y<- x^( -(marginal.post.s2$a+1) ) * exp( logconst-marginal.post.s2$b/x )


plot(x,y,type='l', xlab='Variance of Temperature', ylab='Posterior density')

postscript('chap1_plot3.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
plot(x,y,type='l', xlab='Variance of Temperature', ylab='Posterior density')
graphics.off()


# ---------------------------------------
#
# Plots for the joint posterior of mu and tau
# ---------------------------------------



xrange<-3
x1 <- seq( marginal.post.mu$mu - xrange*marginal.post.mu$sd, marginal.post.mu$mu + xrange*marginal.post.mu$sd, 0.1*marginal.post.mu$sd )

xrange<-3
x2 <- seq( max(c(0, marginal.post.tau$mean - xrange*marginal.post.tau$sd)), marginal.post.tau$mean + xrange*marginal.post.tau$sd, 0.1*marginal.post.tau$sd)

d1<-length(x1)
d2<-length(x2)

mat1<-matrix( x1, nrow=d1,ncol=d2 )
mat2<-matrix( x2, nrow=d1,ncol=d2, byrow=TRUE )

xx1<-as.vector(mat1)
xx2<-as.vector(mat2)
yy <- dnorm( xx1, full.post$mu, sqrt( full.post$c / xx2 )) * dgamma(  xx2, full.post$a, full.post$b )

y<-matrix(yy, nrow=d1,ncol=d2 )


postx1<-x1
postx2<-x2
postymat<-y

contour(postx1,postx2,postymat, xlab='Mean Temperature', ylab='Precision of Temperature')

persp(postx1,postx2,postymat, d=10, theta=30, phi=40)

postscript('chap1_plot4.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par( xaxs='i', yaxs='i', bty='l' , cex=1.5)
contour(postx1,postx2,postymat, xlab='Mean Temperature', ylab='Precision of Temperature')
graphics.off()

postscript('chap1_plot5.ps', width = 10.0, height = 6.0, horizontal=FALSE)
persp(postx1,postx2,postymat, d=10, theta=30, phi=40, xlab='Mean Temperature', ylab='Precision of Temperature', zlab='Joint Posterior Density')
graphics.off()


# ---------------------------------------
#
# Plots for the joint prior of mu and s2
# ---------------------------------------

prior$sd<- prior$c/(prior$b*prior$a)

xrange<-3
x1 <- seq( prior$mu - xrange*prior$sd, prior$mu + xrange*prior$sd, 0.1*prior$sd )

xrange<-3

x2 <- seq( max(c(0, prior$a/prior$b - xrange*sqrt(prior$a)/prior$b)), prior$a/prior$b + xrange*sqrt(prior$a)/prior$b, 0.1*sqrt(prior$a)/prior$b)

d1<-length(x1)
d2<-length(x2)

mat1<-matrix( x1, nrow=d1,ncol=d2 )
mat2<-matrix( x2, nrow=d1,ncol=d2, byrow=TRUE )

xx1<-as.vector(mat1)
xx2<-as.vector(mat2)
yy <- dnorm( xx1, full.post$mu, sqrt( full.post$c / xx2 )) * dgamma(  xx2, full.post$a, full.post$b )

y<-matrix(yy, nrow=d1,ncol=d2 )


priorx1<-x1
priorx2<-x2
priorymat<-y

rangex1<-range(c(priorx1,postx1))
rangex2<-range(c(priorx2,postx2))
rangex3<-range(c(priorymat,postymat))


contour(priorx1,priorx2,priorymat, xlab='Mean Temperature', ylab='Precision of Temperature', xlim=rangex1, ylim=rangex2)

persp(priorx1,priorx2,priorymat, xlab='Mean Temperature', ylab='Precision of Temperature', xlim=rangex1, ylim=rangex2,d=10, zlim=rangex3, theta=30, phi=40)


graphics.off()