Hi all,
I implemented componentwise adaptive Metropolis (Heikki Haario, et
al.) method in R, but it gives different results compared to WinBUGS
simulation rezults.
E. g. suppose you have poisson rate trend like 1.2*1.2^t (where t=1:6
years). I simulated 6 random numbers (1 2 2 1 4 4) and tryint to fit
regression curve a*b^t.
WinBUGS mean for parameter a is around 0.98 and for b is about 1.3.
But adaptive Metropolis (in R) gives posterior mean for parameter a
around 1.5 and for parameter b around 1.19.
But the curve 1.5*1.19^t seems to have poor fit, because all data
poinst (except one is under the curve). And it seems quite disturbing
that estimated curve lies too far from the data points.
Could it be because of different MCMC methods used (some of the
methods in winbugs and adaptime metropolis), or it is just the problem
of my R codding error? (I supply my code below).
I hope, that the problem i encountered is clear and that some of you
might help me.
data=read.csv('failures.csv',header=T,sep=';');
attach(data)
data=data.matrix(data)
#loglikelihood function for poisson regression
loglike=function(a,b){
t=c(1:6)
lamb=a*exp(b*t)
band=data*log(lamb)
suma=sum(band)-sum(lamb)-sum(lfactorial(data))
suma
}
N=100000
Sd=1
eps=0.00000000000001
X=rep(0,n)
X=matrix(data=X,nrow=2,ncol=(n/2),byrow=TRUE)
X[,1]=runif(2,min=0.01,max=0.03)
#variances of parameters
C=rep(0,2)
C[1]<-1; C[2]<-0.1;
T0=100
#non adaptive phase
for(i in 2:T0){
Y=(mvrnorm(1,X[1,i-1],C[1]))
while(Y<0){ Y=(mvrnorm(1,X[1,i-1],C[1]))}
rho=loglike(Y,X[2,i-1])+log(dunif(Y,0,1000))-(loglike(X[1,i-1],X[2,i-1])+log(dunif(X[1,i-1],0,1000)))
X[1,i]=X[1,i-1]+(Y-X[1,i-1])*(log(runif(1))<rho)
Y=(mvrnorm(1,X[2,i-1],C[2]))
while(Y<0){ Y=(mvrnorm(1,X[2,i-1],C[2]))}
rho=loglike(X[1,i],Y)+log(dunif(Y,0,1000))-(loglike(X[1,i],X[2,i-1])+log(dunif(X[2,i-1],0,1000)))
X[2,i]=X[2,i-1]+(Y-X[2,i-1])*(log(runif(1))<rho)
}
#Adaptive phase
mean1=rowMeans(X[1:2,1:(i-2)])
for(i in (T0+1):N/2){
mean1=((i-3)/(i-2))*mean1+X[1:2,(i-2)]/(i-2)
mean2=((i-2)/(i-1))*mean1+X[1:2,(i-1)]/(i-1)
C[1]<-((i-3)/(i-1))*C[1]+Sd*((mean1[1]-mean2[1])^2+((X[1,i-2]-mean2[1])^2)/(i-1))
C[2]<-((i-3)/(i-1))*C[2]+Sd*((mean1[2]-mean2[2])^2+((X[2,i-2]-mean2[2])^2)/(i-1))
Y=(mvrnorm(1,X[1,i-1],C[1]))
while(Y<0){ Y=(mvrnorm(1,X[1,i-1],C[1]))}
rho=loglike(Y,X[2,i-1])+log(dunif(Y,0,1000))-(loglike(X[1,i-1],X[2,i-1])+log(dunif(X[1,i-1],0,1000)))
X[1,i]=X[1,i-1]+(Y-X[1,i-1])*(log(runif(1))<rho)
Y=(mvrnorm(1,X[2,i-1],C[2]))
while(Y<0){ Y=(mvrnorm(1,X[2,i-1],C[2]))}
rho=loglike(X[1,i],Y)+log(dunif(Y,0,1000))-(loglike(X[1,i],X[2,i-1])+log(dunif(X[2,i-1],0,1000)))
X[2,i]=X[2,i-1]+(Y-X[2,i-1])*(log(runif(1))<rho)
}
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
