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