#---------------------------------------------------------------------
#  Example 3.1 

par(mfcol=c(1,2))


theta.vec <- seq(-10,10,length=20001)
theta.mid.int <- theta.vec[-1] - 0.001 / 2


#  1. mu_0 = 0, tau^2_0 = 4, sigma^2 = 1, y = 3


theta.prior   <- diff(pnorm(theta.vec,sd=2))
theta.lik     <- dnorm(3,mean=theta.mid.int,sd=1)
theta.post    <- theta.prior * theta.lik

theta.prior <- theta.prior / sum(theta.prior) * 1000
theta.lik   <- theta.lik   / sum(theta.lik) * 1000
theta.post  <- theta.post  / sum(theta.post) * 1000

plot(theta.mid.int,theta.prior,ylab="Density",xlab="Theta",type="l",xlim=c(-4,6),ylim=c(0,.6),lty=1,col=1)
lines(theta.mid.int,theta.lik,lty=3,col=3)
lines(theta.mid.int,theta.post,lty=4,col=4)

theta.post  <- theta.post  / sum(theta.post) 
sum(theta.post * theta.mid.int)  #  =   (0/4  + 3/1) / ( 1/4 + 1) 
sum(theta.post * theta.mid.int^2) - sum(theta.post * theta.mid.int)^2  #  = 1 / (1/4 + 1/1)


#  2. mu_0 = 0, tau^2_0 = 1, sigma^2 = 1, y = 3

theta.prior   <- diff(pnorm(theta.vec,sd=1))
theta.lik     <- dnorm(3,mean=theta.mid.int,sd=1)
theta.post    <- theta.prior * theta.lik

theta.prior <- theta.prior / sum(theta.prior) * 1000
theta.lik   <- theta.lik   / sum(theta.lik) * 1000
theta.post  <- theta.post  / sum(theta.post) * 1000

plot(theta.mid.int,theta.prior,ylab="Density",xlab="Theta",type="l",xlim=c(-4,6),ylim=c(0,.6),lty=1,col=1)
lines(theta.mid.int,theta.lik,lty=3,col=3)
lines(theta.mid.int,theta.post,lty=4,col=4)

theta.post  <- theta.post  / sum(theta.post) 
sum(theta.post * theta.mid.int)  #  =   (0/1  + 3/1) / ( 1/1 + 1) 
sum(theta.post * theta.mid.int^2) - sum(theta.post * theta.mid.int)^2  #  = 1 / (1/1 + 1/1)



#---------------------------------------------------------------------
#  Example 3.2

 
library(LearnBayes)

?marathontimes

# marathontimes data: 
# Running times in minutes for twenty male runners between
# the ages 20 and 29 who ran the New York Marathon. 

data(marathontimes)
attach(marathontimes)

(s.2   <- var(time))
(n     <- length(time))
(y.bar <- mean(time))

sigma.2 <- s.2 / ( rchisq(1000,n-1) / (n-1))
mu	  <- rnorm(1000,mean=y.bar,sd = sqrt(sigma.2/n)) 	

plot(mu,sigma.2)
points(y.bar,s.2,lty=2,col=2,pch=3,lwd=3)

par(mfcol=c(1,2))
hist(mu)
hist(sigma.2)



#  construct posterior dist of mu 


sigma.2 <- s.2 / ( rchisq(100000,n-1) / (n-1))
mu.smp  <- rnorm(100000,mean=y.bar,sd = sqrt(sigma.2/n)) 	

range(mu.smp)
mu.ser  <- seq(200,350,by = 0.01)
k       <- length(mu.ser)
mu.cdf  <- (rank(c(mu.ser,mu.smp))[1:15001] - (1:15001)) / 10^5

par(mfcol=c(1,1))
plot(mu.ser,mu.cdf,type="l")
abline(h=c(0.025,0.5,0.975),lty=2,col=3)
approx(mu.cdf,mu.ser,xout = c(0.025,0.5,0.975))$y

y.bar 
y.bar + qt(0.975,19) * sqrt(s.2 / 20)
y.bar - qt(0.975,19) * sqrt(s.2 / 20)

mu.diff <- mu.ser[-1] - 0.05
mu.pdf  <- diff(mu.cdf)
plot(mu.diff,mu.pdf,type="l")
sum(mu.diff * mu.pdf)



# Verify that posterior dist of (mu - y.bar) / sqrt(s^2 / n) is t 19 df

sigma.2 <- s.2 / ( rchisq(100000,n-1) / (n-1))
mu	  <- rnorm(100000,mean=y.bar,sd = sqrt(sigma.2/n)) 	

t.smp <- (mu - y.bar) / sqrt(s.2/n)

t.ser <- seq(-8,8,by = 0.01)
r.ser <- rank(c(t.ser,t.smp))[1:1601] - (1:1601)
r.ser <- r.ser / 10^5

par(mfcol=c(1,1))

plot(t.ser,r.ser,type="l")
lines(t.ser,pt(t.ser,df=19),col=2)