# -------------------------------------------------WINBUGS CODE
# MODEL
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.01, 0.01)
}

#DATA
list(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=27)

#FIRST INITIAL VALUES
list(alpha=1,beta=1,tau=1,gamma=0.9)

#SECOND INITIAL VALUES ARE RANDOMLY DRAWN FROM 
#THE PRIOR DISTRIBUTIONS

# -------------------------------------------------READING WINBUGS OUTPUT INTO R
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)
plot(x,y)

#Reading WINBUGS output (2 MCMC chains of length 10000)
M      = 10000
chain1 = matrix(scan("chain1.txt"),4*M,2,byrow=T)
chain2 = matrix(scan("chain2.txt"),4*M,2,byrow=T)
chain1 = matrix(chain1[,2],M,4)
chain2 = matrix(chain2[,2],M,4)

names = c("alpha","beta","gamma","tau")
par(mfrow=c(2,2))
for (i in 1:4){
  ts.plot(chain1[,i],xlab="",ylab="",main="")
  title(names[i])
  lines(chain2[,i],col=2)
}