#---------------------------------------------------------------------
#  Example 2.1 


(0.005 * 0.997) / ( 0.005 * 0.997 + 0.999 * 0.003)




#---------------------------------------------------------------------
#  Example 2.3


theta.vec <- seq(-10,10,length=20001)

plot(theta.vec,theta.vec^2,xlim=c(-3,3),ylim=c(0,3),type="l",xlab="Theta",ylab="Frequentist risk")
lines(theta.vec,.5^2+(1-.5)^2*theta.vec^2,lty=2,col=2)
lines(theta.vec,.8^2+(1-.8)^2*theta.vec^2,lty=3,col=3)
lines(theta.vec,1^2+(1-1)^2*theta.vec^2,lty=4,col=4)
lines(theta.vec,1.2^2+(1-1.2)^2*theta.vec^2,lty=5,col=5)
lines(theta.vec,1.3^2+(1-1.3)^2*theta.vec^2,lty=6,col=6)


 

# Numerical computation of posterior for N(0,1) prior 

theta.vec[1:4]
diff(theta.vec[1:4])

theta.mid.int <- theta.vec[-1] + 0.001 / 2
theta.prior   <- diff(pnorm(theta.vec))

#  likelihood and posterior for y = 1

y <- 1

theta.lik  <- dnorm(y,mean=theta.mid.int,sd=1)
theta.post <- theta.prior * theta.lik

plot(theta.mid.int,theta.prior,ylab="Density",xlab="Theta",type="l",xlim=c(-4,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.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)

#  likelihood and posterior for y = 2

y <- 2

theta.lik  <- dnorm(y,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)

lines(theta.mid.int,dnorm(theta.mid.int,mean=1,sd=sqrt(1/2)),col=4)


#  Compute numerical mean and variance

sum(theta.post * theta.mid.int)


# First posterior must be turned into a discrete distribution!

theta.post <- theta.post / sum(theta.post)

(mn.theta <- sum(theta.post * theta.mid.int))
sum(theta.post * theta.mid.int^2) - mn.theta^2



# Derive 0.95 credible interval

(ci.95 <- qnorm(c(0.025,0.975),mean=1,sd=sqrt(1/2)))

arrows(ci.95[1],0,ci.95[2],0,angle=90,code=3,length=0.05,col=4)
lines(theta.mid.int,dnorm(theta.mid.int,mean=1,sd=sqrt(1/2)),col=4)


plot(theta.vec,cumsum(c(0,theta.post))/1000,type="l",ylab="Cumulative distribution",xlab="Theta")

plot(theta.vec,cumsum(c(0,theta.post))/1000,type="l",lty=2,col=4,xlim=c(-1,3),ylab="Cumulative distribution",xlab="Theta")
lines(theta.vec,pnorm(theta.vec,mean=1,sd=sqrt(1/2)),lty=1,col=4)

abline(h=0.025,col=3) ; abline(h=0.975,col=3)

abline(v=ci.95[1],col=4,lty=2) ; abline(v=ci.95[2],col=4,lty=2)
arrows(ci.95[1],0,ci.95[2],0,angle=90,code=3,length=0.05,col=4)

approx(theta.vec,cumsum(c(0,theta.post))/1000,xout=c(0))
pnorm(0,mean=1,sd=sqrt(1/2))

approx(cumsum(c(0,theta.post))/1000,theta.vec,xout=c(0.025,0.975))$y
ci.95