A colleague who is very good at C++ has coded a Metropolis-Hastings
algorithm to estimate the posteriors of a very complex spatial statistical
model.
He has tuned his code so that it has an acceptance rate of around 30%, but
the
results we are obtaining do not match those of WinBUGS. I have pasted
the relevant WinBUGS code below, and any comments would be welcome.
Based on our WinBUGS results, we are convinced that WinBUGS is giving
us realistic posterior estimates since they correspond very well to some
kriging results obtained using Splus.
I would be especially interested in learning if a Slice Sampler or some
other
estimation method would be more preferred. Any suggestions as to how to
duplicate
this code with a C++ MCMC method would be GREATLY appreciated!
Thanks in advance!
David Paul, Ph.D.
Battelle Memorial Institute
505 King Avenue
Columbus, OH 43201
614.424.3176
model
{
#First, we define the errors in the SAR model structure for the
#data. The term 'SAR' is used by Ripley in his Spatial Statistics textbook
#and denotes a 'simultaneous autoregressive' model. The implementation
#here is slightly different from that of Ripley and Cressie. The
autoregressive
#terms are coded very similarly to that found in Peter Congdon's new book
#on Bayesian Statistical Modelling. The application in question deals
#with a regular lattice, which makes a SAR model attractive.
for(i in 1:Num.Total.LatticeGrids)
{
e[i] ~ dnorm(e.bar[i],tau.Error)
#The following line specifies that the mean of the autoregressive "errors"
of all neighboring
#sites should be multiplied by rho. This quanity forms the mean of the
distribution of the
#autoregressive "error" at the current Lattice Grid Cell.
e.bar[i] <- rho*mean( e.neigh[
(cumulative.neighbor.number[i]+1):cumulative.neighbor.number[i+1] ] )
}
#This line causes WinBUGS to associate the numbering of the spatial
autoregressive
#"errors" with the numbering of the Grid Cells in the Regular Lattice.
for(m in 1:Num.Total.Neighbors)
{
e.neigh[m] <- e[neighbor[m]]
}
#Now we define the likelihood of the data. We are going to impose a SAR
#structure on the Predictive errors.
for(i in 1:Num.Total.Grids)
{
#X[i,1] represents the latitude of an observation and X[i,2] represents the
#longitude of an observation, so mu.Trend is a simple trend surface with a
spatial
#autoregressive term added to the end.
mu.Trend[i]<-Beta[1] + Beta[2]*X[i,1] + Beta[3]*X[i,2] + e[i]
#There are two sets of data that we are using to obtain estimates of the
posterior
#spatial surface. The variable "cumulative.observations" is an index for
the total
#number of observations within a Grid Cell on the Regular Lattice. It is
assumed
#that one of the data sets is biased upwards from the truth, so that
mu.Trend needs
#to be appended with a constant term, a. The function equals( ) is designed
to decide
#if a particular observation belongs to the data set needing the adjustment.
#It is also assumed that the precision associated with each data set is
different, so
#that tau changes as well.
for(j in (cumulative.observations[i]+1):cumulative.observations[i+1])
{
logic[j] <- a*equals(indicator[j],2)
mu.Observation[j] <- logic[j] + mu.Trend[i]
Y.Combined[j] ~ dnorm(mu.Observation[j],tau[indicator[j]])
}
}
#PRIORS, not including the autoregressive terms - they are defined above
tau.Error ~ dgamma(0.0001,0.0001)
for(i in 1:2)
{
tau[i] ~ dgamma(0.0001,0.0001)
}
Beta[1:3] ~ dmnorm(Beta.Mean[] , Beta.R[,])
Beta.R[1:3,1:3] ~ dwish(Beta.Omega[,] , 3)
a ~ dnorm(0.0 , 0.001)
rho ~ dunif(0,1)
}
-------------------------------------------------------------------
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
|
