#########################################################################
#
# y[i] : velocity of an enzymatic reacton (in counts/min/min) 
#        as a function of substrate concentration
# x[i] : (in ppm) where the enzyme has been treated with Puromycin
# 
# Referencia: Carlin and Louis (1998) - Bayes and Empirical 
# Bayes Methods for Data Analysis - Chapman & Hall/CRC
# paginas: 233-234.
#
# Hedibert F. Lopes
# www.hedibert.org
# February 2nd 2021
#
#########################################################################

mu       = function(th){gamma+alpha*x/(th+x)}
logprior = function(th){dnorm(th,theta0,Ctheta,log=TRUE)}
loglike  = function(th){sum(dnorm(y,mu(th),sqrt(tau2),log=TRUE))}
logpost  = function(th){dnorm(th,theta0,Ctheta,log=TRUE)+loglike(th)}

# Data
n        = 12
x        = c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
y        = c(76,47,97,107,123,139,159,152,191,201,207,200)
plot(x,y,xlab="Substrate concentration (in ppm)",
     ylab="Velocity of an enzymatic reacton (in counts/min/min)",
     pch=16,cex=1.5)

# Prior, likelihood and posterior
gamma    = 50
alpha    = 170
tau2     = 126
theta0   = 0.0
Ctheta   = 100
N        = 1000
thetas   = seq(0.08,0.20,length=N)
logposts = rep(0,N)
for (i in 1:N)
  logposts[i] = logpost(thetas[i])
posts   = exp(logposts)
posts   = posts/max(posts)
thetas1 = seq(-5,5,length=N)
priors  = dnorm(thetas1,theta0,sqrt(Ctheta))


par(mfrow=c(1,2))
plot(thetas1,dnorm(thetas1,theta0,sqrt(Ctheta)),type="l",ylab="Prior",
     xlab="theta")
abline(v=0.08,col=2)
abline(v=0.2,col=2)
plot(thetas,posts,type="l",xlim=c(0,0.3),xlab="theta",
     ylab="Posterior (up to a normalizing constant")
lines(thetas1,priors/max(priors),col=2)
     


# SIR-based posterior inference for theta
set.seed(655432)
M    = 10000
th.d = runif(M,0.0,0.3)
w    = rep(0,M)
for (i in 1:M)
  w[i] = logpost(th.d[i])
w = exp(w-max(w))
w = w/sum(w)
th.sir = sample(th.d,replace=TRUE,size=M,prob=w)

par(mfrow=c(1,3))
hist(th.d,main="Draws from proposal U(0,0.3)",xlab="theta",prob=TRUE)
hist(th.sir,main="SIR draws from posterior",xlab="theta",prob=TRUE)
plot(thetas,exp(logposts),type="l")
title("Posterior of theta\n up to a normalizing constant")

# Summaries from the posterior of theta
mean(th.sir)
sqrt(var(th.sir))
quantile(th.sir,c(0.005,0.995))


# Random-Walk-Metropolis-Hastings (RW-MH)-based 
# posterior inference for theta
set.seed(192796)
theta  = 2
Vtheta = (0.025)^2
burnin = 1000
Lag    = 1
M      = 10000
niter  = burnin+Lag*M
th.rw  = rep(0,niter)
for (iter in 1:niter){
  theta1 = rnorm(1,theta,sqrt(Vtheta))
  acceptance = exp(logpost(theta1)-logpost(theta))
  if (runif(1) < acceptance){
    theta = theta1
  }
  th.rw[iter] = theta
}
th.rw = th.rw[seq(burnin+1,niter,by=Lag)]

# RW-MH - trace plot and ACF for draws
par(mfrow=c(2,1))
ts.plot(th.rw,type="l",xlab="Iterations",ylab="theta")
acf(th.rw,main="")

# Posterior densities via RW-MH
par(mfrow=c(1,2))
hist(th.rw,prob=TRUE,xlab="theta",main="RW-MH draws",ylab="Posterior density")
plot(thetas,exp(logposts),type="l",xlab="theta",ylab="Posterior density up to a normalizing constant")




# Comparing SIR and RW-MH
par(mfrow=c(2,2))
ts.plot(th.rw,xlab="Iterations",ylab="theta",main="RW-MH draws")
acf(th.rw,main="")
hist(th.rw,main="RW-MH draws",xlab="theta",prob=TRUE,ylab="Posterior density")
hist(th.sir,main="SIR draws",xlab="theta",prob=TRUE,ylab="Posterior density")

par(mfrow=c(1,1))
plot(density(th.sir),xlab="theta",main="",lwd=3,ylab="Posterior density")
lines(density(th.rw),col=2,lwd=3)
legend("topright",legend=c("SIR","RW-MH"),col=1:2,lwd=3,bty="n")

tab = rbind(c(mean(th.sir),sqrt(var(th.sir)),quantile(th.sir,c(0.005,0.995))),
c(mean(th.rw),sqrt(var(th.rw)),quantile(th.rw,c(0.005,0.995))))

rownames(tab) = c("SIR","RW-MH")
colnames(tab) = c("Mean","StDev","0.5%","99.5%")
round(tab,4)