#source("c:\\David\\Enseignements\\Cours ICASA\\MHresponse.txt")
####Estimation of three parameters of a response curve using a Metropolis-Hastings algorithm#### 
#Makowski D, Lavielle M. 2006.Using SAEM to estimate parameters of models of response to applied fertilizer. JABES 11, 45-60.
###Seed used for random sampling### 
set.seed(1)
###Number of simulations### 
N<-50000
###Prior distribution parameters### 
#Maximum yield 
mu.1<-9.18 
tau.1<-1.16 
#Crop requirement 
mu.2<-0.026 
tau.2<-0.0065 
#Optimal dose 
mu.3<-123.85 
tau.3<-46.70
###Data### 
Y<-c(2.50, 5.01, 7.45, 7.51) 
#Y<-c(7.51, 9.6, 10.11, 10.10) 
X<-c(0, 50, 100, 200) 
sigma<-0.3*100
###Initialisation of the parameter chains### 
Theta.1<-1:N 
Theta.1[1]<-mu.1 
Theta.2<-1:N 
Theta.2[1]<-mu.2 
Theta.3<-1:N 
Theta.3[1]<-mu.3
###Computation of the likelihood for the first set of parameter values### 
Likeli<-matrix(nrow=N,ncol=length(Y)) 
LikeliProd<-1:N 
Reject<-0
for(j in 1:length(Y)) 
{ 
Ymax.i<-Theta.1[1] 
B.i<-Theta.2[1] 
Xmax.i<-Theta.3[1] 
Model.i.j<-min(Ymax.i, Ymax.i+B.i*(X[j]-Xmax.i)) 
Likeli[1,j]<-max(dnorm(Y[j],Model.i.j,sigma), 10^(-100)) 
}
LikeliProd[1]<-max(prod(Likeli[1,]), 10^(-100)) 
###Tuning parameters### 
Lambda.1<-0.3 
Lambda.2<-0.3 
Lambda.3<-0.3
###Generation of candidate parameter values and test###
for (i in 2:N) {
Theta.1c<-rnorm(1,Theta.1[i-1],Lambda.1*tau.1) 
Theta.2c<-rnorm(1,Theta.2[i-1],Lambda.2*tau.2) 
Theta.3c<-rnorm(1,Theta.3[i-1],Lambda.3*tau.3)
for(j in 1:length(Y)) { 
Ymax.i<-Theta.1c 
B.i<-Theta.2c 
Xmax.i<-Theta.3c 
Model.i.j<-min(Ymax.i, Ymax.i+B.i*(X[j]-Xmax.i)) 
Likeli[i,j]<-max(dnorm(Y[j],Model.i.j,sigma),10^(-100)) 
} 
LikeliProd[i]<-max(prod(Likeli[i,]), 10^(-100)) 

Test.i<-log(LikeliProd[i]) + log(dnorm(Theta.1c,mu.1,tau.1)) + log(dnorm(Theta.2c,mu.2,tau.2)) + log(dnorm(Theta.3c,mu.3,tau.3))- 
log(LikeliProd[i-1])-log(dnorm(Theta.1[i-1],mu.1,tau.1))-log(dnorm(Theta.2[i-1],mu.2,tau.2))-log(dnorm(Theta.3[i-1],mu.3,tau.3))
u<-runif(1,min=0,max=1) 
if (Test.i>log(u)) { 
Theta.1[i]<-Theta.1c 
Theta.2[i]<-Theta.2c 
Theta.3[i]<-Theta.3c 
}
if (Test.i<=log(u)) { 
Theta.1[i]<-Theta.1[i-1] 
Theta.2[i]<-Theta.2[i-1] 
Theta.3[i]<-Theta.3[i-1] 
Reject<-Reject+1 
} 
} 
###Chains of parameter values 
par(mfrow=c(2,2)) 
plot(1:N, Theta.1, xlab="Iteration", ylab="Theta 1", type="l", cex=2) 
plot(1:N, Theta.2, xlab="Iteration", ylab="Theta 2", type="l", cex=2) 
plot(1:N, Theta.3, xlab="Iteration", ylab="Theta 3", type="l", cex=2)
readline()
###Posterior distribution### 
par(mfrow=c(2,2)) 
hist(Theta.1[1:(N/2)],labels=F,xlab="Theta.1",main="Posterior distribution of Theta.1",cex=2) 
hist(Theta.2[1:(N/2)],labels=F,xlab="Theta.2",main="Posterior distribution of Theta.2",cex=2) 
hist(Theta.3[1:(N/2)],labels=F,xlab="Theta.3",main="Posterior distribution of Theta.3",cex=2)
print("Posterior mean:") 
print(mean(Theta.1[1:(N/2)])) 
print(mean(Theta.2[1:(N/2)])) 
print(mean(Theta.3[1:(N/2)])) 
print("Posterior variance:") 
print(var(Theta.1[1:(N/2)])) 
print(var(Theta.2[1:(N/2)])) 
print(var(Theta.3[1:(N/2)])) 
print("Acceptance rate (%):") 
print(100*(N-Reject)/N) 

###Plot of average and estimated response curves### 
Dose<-1:200 
Average.Response<-mu.1+mu.2*(Dose-mu.3) 
Average.Response[Average.Response>mu.1]<-mu.1 
Estimated.Response<-mean(Theta.1[1:(N/2)])+mean(Theta.2[1:(N/2)])*(Dose-mean(Theta.3[1:(N/2)])) 
Estimated.Response[Estimated.Response>mean(Theta.1[1:(N/2)])]<-mean(Theta.1[1:(N/2)])
plot(Dose, Average.Response, xlab="N dose (kg/ha)", ylab="Yield (t/ha)", type="l", ylim=c(2,12),cex=2) 
lines(Dose, Estimated.Response, lty=4, lwd=2) 
points(X,Y,pch=20,cex=2)