#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)