##################################################
# Example 3.6 from Gamerman and Lopes (2006)
##################################################

# Computing the posterior mean and posterior variance of theta
set.seed(12154)
M = 10000
theta = runif(M)

# Posterior mean via MC integration
numerator = mean((2+theta)^(125)*(1-theta)^(38)*theta^(35))
denominator = mean((2+theta)^(125)*(1-theta)^(38)*theta^(34))
Etheta.mc = numerator/denominator


# Posterior variance via MC integration
numerator1 = mean((theta-Etheta.mc)^2*(2+theta)^(125)*(1-theta)^(38)*theta^(34))
Vtheta.mc = numerator1/denominator

# Approximated 95% posterior interval (approximated because this is not the correct computation)
c(Etheta.mc-2*sqrt(Vtheta.mc),Etheta.mc+2*sqrt(Vtheta.mc))


# Obtaining draws from the posterior p(theta) via SIR
w = (2+theta)^(125)*(1-theta)^(38)*theta^(34)
plot(theta,w)

theta1 = sample(theta,replace=TRUE,prob=w,size=5000)

hist(theta1)

# Posterior mean and variance of theta via SIR
Etheta.sir = mean(theta1)
Vtheta.sir = var(theta1)


# Comparing MC integration and SIR when approximating the posterior mean of theta
# Rao-blackwell theorem in action
set.seed(12154)
M = 10000
replication = 1000
estimators = matrix(0,replication,2)
for (j in 1:replication){
  theta = runif(M)  
  numerator = mean((2+theta)^(125)*(1-theta)^(38)*theta^(35))  
  denominator = mean((2+theta)^(125)*(1-theta)^(38)*theta^(34))
  Etheta.mc = numerator/denominator
  w = (2+theta)^(125)*(1-theta)^(38)*theta^(34)
  theta1 = sample(theta,replace=TRUE,prob=w,size=5000)
  Etheta.sir = mean(theta1)
  estimators[j,] = c(Etheta.mc,Etheta.sir)
}
boxplot(estimators[,1],estimators[,2],names=c("MC integration","SIR"),ylab="Posterior mean",main="")
title("MC approximations to the posterior mean of theta")










##################################################
# Example 3.7 from Gamerman and Lopes (2006)
##################################################
quantile025 = function(x){quantile(x,0.025)}
quantile975 = function(x){quantile(x,0.975)}
post = function(th){sum((y-th[1]*(1-exp(-x*th[2])))^2)^(-2)}


# Data
# ----
x = c(1,2,3,4,5,7)
y = c(8.3,10.3,19,16,15.6,19.8)
#y = c(8.3,10.3,19,16,15.6,49.8)

# Contours of the posterior distribution
# --------------------------------------
M   =  200
th1 =  seq(-20,50,length=M)
th2 =  seq(-2,10,length=M)
#th1 =  seq(-500,600,length=M)
#th2 =  seq(-0.5,0.5,length=M)
f   =  matrix(0,M,M)
for (i in 1:M)
  for (j in 1:M)
    f[i,j] = post(c(th1[i],th2[j]))
f = f/sum(f)

par(mfrow=c(1,1))
contour(th1,th2,f,nlevels=50,drawlabels=FALSE)


# Sampling importance resampling (SIR) algorithm
# ----------------------------------------------
set.seed(9023749)
M       = 50000
M1 = M/10
ths     = matrix(0,M,2)
ths[,1] = -20+70*runif(M)
ths[,2] = -2+12*runif(M)
#ths[,1] = -500+1100*runif(M)
#ths[,2] = -0.5+1*runif(M)
w = NULL
for (i in 1:M)
  w = c(w,post(c(ths[i,1],ths[i,2])))
w    = w/sum(w)
ind  = sample(1:M,size=M1,replace=T,prob=w)
ths1 = ths[ind,]

# Posterior summaries
# -------------------
matrix(c(apply(ths1,2,mean),
sqrt(apply(ths1,2,var)),
apply(ths1,2,quantile025),
apply(ths1,2,quantile975)),2,4)

par(mfrow=c(1,1))
plot(ths1,xlab=expression(theta[1]),ylab=expression(theta[2]),xlim=range(th1),ylim=range(th2),axes=F)
axis(1)
axis(2)
contour(th1,th2,f,drawlabels=F,nlevels=50,axes=F,add=T,col=2)


# Posterior predictive
# -----------------------
N  = 20
x1 = seq(min(x),max(x),length=N)
f1 = matrix(0,M1,N)
for (i in 1:N)
  f1[,i] = ths1[,1]*(1-exp(-x1[i]*ths1[,2]))
q050 = apply(f1,2,median)
q025 = apply(f1,2,quantile025)
q975 = apply(f1,2,quantile975)


par(mfrow=c(1,1))
plot(x,y,ylim=range(q025,q975),xlab="time (days)",ylab="BOD (mg/liter)",main="",axes=F)
axis(1)
axis(2)
lines(x1,q025,lty=2)
lines(x1,q050)
lines(x1,q975,lty=2)