#######################################################################################
#
# Carlin and Gelfand (1991) present a nonconjugate Bayesian analysis of the 
# following data set from Ratkowsky (1983):
#
#	Dugong	       1     2	   3	 4    5	  ....    26    27
#	__________________________________________________________
#	Age (X)	     1.0   1.5	 1.5   1.5   2.5  ....	29.0  31.5
#	Length (Y)  1.80  1.85	1.87  1.77  2.02  ....	2.27  2.57
#
# The data are length and age measurements for 27 captured dugongs (sea cows). 
# Carlin and Gelfand (1991) model this data using a nonlinear growth curve with no 
# inflection point and an asymptote as Xi tends to infinity:
#
#	Yi  ~  Normal(mi, t),        	i = 1,...,27
#	
#	mi  =  a - b*g^Xi    	a, b > 1; 0 < g < 1
#	
# Standard noninformative priors are adopted for a, b and t, and a uniform prior 
# on (0,1) is assumed for g. However, this specification leads to a non conjugate 
# full conditional distribution for g which is also non log-concave. 
#
#######################################################################################

# BEFORE ANYTHING ELSE, SAVE THE FOLLOWING MODEL IN FILE "nonlinearmodel.bug"

# BEGIN -- BEGIN -- BEGIN -- BEGIN -- BEGIN -- BEGIN -- BEGIN -- BEGIN -- BEGIN 
# -------------------------------------------------------------------------------
model{
  for(i in 1:N){
    y[i]   ~ dnorm(mu[i], tau)
    mu[i] <- alpha-beta*pow(gamma,x[i])			
  }
  alpha ~ dnorm(0.0, 1.0E-6)
  beta ~ dnorm(0.0, 1.0E-6)
  gamma ~ dunif(0.0, 1.0)
  tau ~ dgamma(0.001, 0.001)
}
# -------------------------------------------------------------------------------
# END -- END -- END -- END -- END -- END -- END -- END -- END -- END -- END


quant025=function(x){quantile(x,0.025)}
quant975=function(x){quantile(x,0.975)}

x = c( 1.0,  1.5,  1.5,  1.5, 2.5,   4.0,  5.0,  5.0,  7.0, 8.0,  8.5,  9.0,  9.5, 9.5,  
      10.0, 12.0, 12.0, 13.0, 13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5)

y = c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47, 2.19, 2.26, 2.40, 2.39, 2.41, 
      2.50, 2.32, 2.32, 2.43, 2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57)

N       <- length(x)
data    <- list("x","y","N")
inits   <- function(){
  list(alpha = 1, beta = 1, tau = 1, gamma = 0.9)
}

nonlinear.sim = bugs(data,inits,model.file="nonlinearmodel.bug",
                     parameters=c("alpha","beta","tau","gamma"),n.chains=1,
                     n.iter=10000,n.burnin=5000,n.thin=1,
                     bugs.directory="c:/Program Files/WinBUGS14/",codaPkg=FALSE)
print(nonlinear.sim)
apply(nonlinear.sim$sims.matrix,2,summary)

alpha = nonlinear.sim$sims.matrix[,1]
beta  = nonlinear.sim$sims.matrix[,2]
tau   = nonlinear.sim$sims.matrix[,3]
gamma = nonlinear.sim$sims.matrix[,4]


par(mfrow=c(1,2))
xs = seq(min(x),max(x),length=100)
f = NULL
for (i in 1:100)
  f = cbind(f,alpha-beta*gamma^xs[i])
fmean   = apply(f,2,mean)
fmedian = apply(f,2,median)
f025    = apply(f,2,quant025)
f975    = apply(f,2,quant975)
L       = min(c(y,fmean,fmedian,f025,f975))
U       = max(c(y,fmean,fmedian,f025,f975))
plot(x,y,xlab="AGE",ylab="LENGTH",axes=F,ylim=c(1.5,3.0))
axis(1)
axis(2)
lines(xs,fmean,col=2)
lines(xs,fmedian,col=4)
lines(xs,f025,col=3)
lines(xs,f975,col=3)

ys = seq(1.5,3.0,length=50)
xs = seq(min(x),max(x),length=50)
f  = matrix(0,50,50)
for (i in 1:50)
  for (j in 1:50)
    f[i,j] = mean(dnorm(ys[j],alpha-beta*gamma^xs[i],sqrt(1/tau)))
plot(x,y,xlab="AGE of new dugongs",ylab="LENGTH of new dugongs",axes=F,xlim=range(xs),ylim=range(ys))
axis(1)
axis(2)
contour(xs,ys,f,drawlabels=F,add=T,col=2)

par(mfrow=c(1,1))
persp(xs,ys,f,theta=0,phi=90,expand=0.5,col="lightblue",
      ltheta=120,shade=0.75,ticktype="detailed",
      xlab="AGE of new dugongs",ylab="LENGTH of new dugongs")