Tuesday, March 04, 2008

An Example for Just Another Gibbs Sampler (JAGS)

JAGS is a cross platform BUGS. Because it runs in a command mode, it is somewhat faster than other BUGS. However, I found it is a bit tacky to use it as a Windows/PC user. I found the manual very useful but some of the language is a bit uncommon to regular users. For instance, compile [,nthins( < n > )]. Do we need brackets here or do we need <> here? It turns out that we need neither.

So I decide to write this entry to take a note for myself and to help a user like me who is not sophisticated enough to understand the manual completely.

First, download JAGS and install it. For Windows users, the ready-to-install binary is very easy to install. For others, I have learned that it is not an easy task to install JAGS properly.

Second, an example is always helpful to understand how JAGS works! Here is the model (model.txt), the data (jagsdata.txt), the initial values (jagsinits1.txt, jagsinits2.txt, jagsinits3.txt) and the script (jagsscript.txt).

For JAGS, the way to input data and input inits is different. Don't use "=". Use "<-". And each parameters as well as variables are a single vector or a matrix. Take a look at the example files for data and inits, you will understand what I mean.

So once you save all the example file in the same directory, let's say "c:/". So we have to change working directory first (all the codes below are in jagsscript.txt):

cd "c:/"

Then we read in the model and the data, compile the model, read in inits, initialize the model, update the model (as the burin), monitor parameters, update the model again, and save code statistics.

model in "model.txt"
data in "jagsdata.txt"
compile, nchains(3)
inits in "jagsinits1.txt"
inits in "jagsinits2.txt"
inits in "jagsinits3.txt"
initialize
update 1000
monitor mu, thin(5)
monitor sigma, thin(5)
monitor theta, thin(5)
monitor deviance, thin(5)
monitor, type(pD)
update 5000
coda *
monitors to "DEV.txt"
monitors to "pD.txt", type(pD)

You can do the above code in the command mode as:

jags jagsscript.txt

Or you can run the script in the JAGS terminal as:

run jagsscript.txt

Then you can check to see if there are "CODAchain1.txt", "CODAchain2.txt", "CODAchain3.txt" and "CODAindex.txt" under your working directory. They are 3 coda chains values and 1 index file.

To read in these coda files as the BUGS object, please download this file jags2bugs.R and load it into R by source("jags2bugs.R").

After that,

library(R2WinBUGS)
jags.fit <- jags2bugs("c:/", para=c("mu", "sigma", "theta"))
print(jags.fit)
plot(jags.fit)

For some reason, I was not able to get the pD value right here. But in any case, jags2bugs will calculate deviance and pD for you. So we don't need the pD and deviance parts here acutally. The script can be simplified as:

model in "model.txt"
data in "jagsdata.txt"
compile, nchains(3)
inits in "jagsinits1.txt"
inits in "jagsinits2.txt"
inits in "jagsinits3.txt"
initialize
update 1000
monitor mu, thin(5)
monitor sigma, thin(5)
monitor theta, thin(5)
update 5000
coda *

3 comments:

Ine said...

Thank you!

Anonymous said...

Thank you very much

Margarita Rincón Hidalgo said...

Thank you!!, I just want to say that coda archives could be read within R with codamenu() from coda package and it has a lot of tools to analyse convergence.