# program to implement Adaptive Metropolis algorithm of Haario et al. 
#
# Reference:
# An adaptive Metropolis algorithm, by H. Haario, E. Saksman and J. Tamminen,
# Bernoulli 7(2), 2001, 223-242
#
AdaptMH=function(lh,par,npar,C0,nsim,n1,n2,n3,eps,scal){
#
# lh: neg log likelihood function to be called
# par: initial set of parameter estimates
# npar: dimension of par
# C0: npar x npar matrix, initial guess of covariance matrix
# nsim: total number of simulations
# n1: initial run with given C0
# n2: interval between subsequent updates
# n3: interval between updates of output parameter matrix
# scal: scaling parameter, recommended value is 2.4 but user is allowed to change
#   sd parameter is scal*scal/npar, see Haario 2001, Gelman et al 1996
# parout: output parameter matrix of dimension n4 x npar, where n4=floor(nsim/n3)
#
n4=floor(nsim/n3)
parout=matrix(nrow=n4,ncol=npar)
C=C0
sum1=par
sum2=par %*% t(par)
sd=scal*scal/npar
ID=matrix(0,nrow=npar,ncol=npar)
diag(ID)=rep(1,npar)
parold=par
lhold=lh(par)
acc0=0
for(isim in 1:nsim){
if(isim==n1){
# update C and xbar
xbar=sum1/(n1+1)
C=(sum2-(n1+1)*xbar%*%t(xbar))/n1
sum1=rep(0,npar)
sum2=matrix(0,nrow=npar,ncol=npar)
}
if(isim>n1&(n2*floor((isim-n1)/n2)==(isim-n1))){
# update C and xbar
sum2=sum2+(isim-n2+1)*xbar%*%t(xbar)
xbar=((isim-n2+1)*xbar+sum1)/(isim+1)
C=(isim-n2)*C/isim+(sum2-(isim+1)*xbar%*%t(xbar))/isim
sum1=rep(0,npar)
sum2=matrix(0,nrow=npar,ncol=npar)
}
# generate new trial value
parnew=as.vector(parold+t(chol(sd*C+sd*eps*ID))%*%rnorm(npar))
lu=log(runif(1))
lhnew=lh(parnew)
if(lhnew-lhold+lu<0){
parold=parnew
lhold=lhnew
acc0=acc0+1
}
iout=floor(isim/n3)
if(n3*iout==isim){
# add current parameter vertor to output matrix
parout[iout,]=parold
}
sum1=sum1+parold
sum2=sum2+parold%*%t(parold)
}
return(value=list(parout=parout,acc0=acc0))
}