export LINBUGS_HOME="/c/tduong/home/bin/openbugs" export LINBUGS_TEMP="/c/tduong/home/temp"where LINBUGS_HOME is the directory where you want to install LinBugs and LINBUGS_TEMP is the directory where LinBugs writes its temporary files.
http://mathstat.helsinki.fi/openbugsUnpack this into $LINBUGS_HOME. After unpacking, check that there is a file called cbugs. If not, download it from
http://www.gatsby.ucl.ac.uk/~iam23/compnotes/openbugsIf there is already a file called linbugs in $LINBUGS_HOME, copy it to another file e.g. linbugs.orig. Create a file called linbugs containing the following commands
#!/bin/bash export LD_ASSUME_KERNEL=2.4.1 cd $LINBUGS_HOME if [ \! -e $LINBUGS_TEMP ] ; then mkdir $LINBUGS_TEMP fi if [ -e bugs.so ] ; then ./cbugs "$LINBUGS_HOME" "$LINBUGS_TEMP" "/bugs.so" else ./cbugs "$LINBUGS_HOME" "$LINBUGS_TEMP" "/brugs.so" fi
./linbugsSomething like
OpenBUGS ClassicBUGS release 2.1.1 type 'modelQuit()' to quit Bugs>should appear in the command line. Quit from LinBugs by typing modelQuit().
######################################################################### ## Simple R interface for LinBugs ######################################################################### ########################################################################## ## Execute LinBugs from inside R session ## ## Parameters ## file - file of LinBUGS commands ## ## Returns ## Invisible ########################################################################## rlinbugs.exec <- function(file) { system(paste("bash $LINBUGS_HOME/linbugs <", file), ignore.stderr=TRUE) invisible() } ########################################################################## ## Reads CODA output files into R objects ## ## Parameters ## stem - stem for CODA output files ## chain.num - chain number ## ## Returns ## Matrix with 3 columns ## first column is "iter" = iteration number ## subsequent columns are the values of the each iteration of parameters ## in MCMC ######################################################################### rlinbugs.read <- function(stem, chain.num=1) { varnames <- read.table(paste(stem,"CODAindex.txt", sep="")) n.var <- nrow(varnames) n.iter <- varnames[1,3] output <- matrix(0, nrow=n.iter, ncol=(n.var+1)) for (i in 1:n.var) { if (i==1) { output.temp <- read.table(paste(stem,"CODAchain", chain.num, ".txt", sep=""), nrows=n.iter) output[,1:2] <- as.matrix(output.temp) } else if (i>1) { output.temp <- read.table(paste(stem,"CODAchain", chain.num, ".txt", sep=""), skip=(i-1)*n.iter, nrows=n.iter) output[,i+1] <- as.matrix(output.temp)[,2] } } colnames(output) <- c("iter", as.character(varnames[,1])) return(output) }Save these into rlinbugs.R in your working directory. In my case it is, /c/tduong/home/research/bugs.
price = a + b*age
In your working directory, create the following files. Unfortunately I can't quite force LinBUGS to read in files from the working directory (its default is to use $LINBUGS_HOME), so you have to specify the full address of all input/ouput files.
File name | Commands |
---|---|
carsbugs.txt |
modelCheck("/c/tduong/home/research/bugs/carsmodel.txt") modelData("/c/tduong/home/research/bugs/carsdata.txt") modelCompile() modelInits("/c/tduong/home/research/bugs/carsinits.txt") modelUpdate(10000) summarySet("a") summarySet("b") samplesSet("a") samplesSet("b") modelUpdate(10000) summaryStats("*") samplesStats("*") samplesCoda("*","/c/tduong/home/research/bugs/output") modelQuit() |
carsmodel.txt |
model { for( i in 1 : N ) { mu[i] <- a+b*age[i] price[i] ~ dnorm(mu[i], tau) } tau ~ dgamma(0.001, 0.001) sigma <- 1 / sqrt(tau) a ~ dnorm(0, 1.0E-12) b ~ dnorm(0, 1.0E-12) } |
carsinits.txt |
list( a=0, b=0, tau=1) |
carsdata.txt |
list(N=39, age = c(13, 14, 14,12, 9, 15, 10, 14, 9, 14, 13, 12, 9, 10, 15, 11, 15, 11, 7, 13, 13, 10, 9, 6, 11, 15, 13, 10, 9, 9, 15, 14, 14, 10, 14, 11, 13, 14, 10), price = c(2950, 2300, 3900, 2800, 5000, 2999, 3950, 2995, 4500, 2800, 1990, 3500, 5100, 3900, 2900, 4950, 2000, 3400, 8999, 4000, 2950, 3250, 3950, 4600, 4500, 1600, 3900, 4200, 6500, 3500, 2999, 2600, 3250, 2500, 2400, 3990, 4600, 450, 4700)) |
In the samplesCoda command in carsbugs.txt, the second argument is known as the stem, and it creates output files of the form outputCODAchainN.txt where N is the chain number, and outputCODAindex.txt, in the directory /c/tduong/home/research/bugs
Note that graphical functions are not yet (as of 31 Dec 2005) available in LinBugs.
source("rlinbugs.R") rlinbugs.exec(file="carsbugs.txt") cars <- rlinbugs.read(stem="/c/tduong/home/research/bugs/output", chain.num=1)The argument to stem in rlinbugs.read should be exactly the same as the stem in the samplesCoda.
The text output in the R command line is something like
OpenBUGS ClassicBUGS release 2.1.1 type 'modelQuit()' to quit Bugs> model is syntactically correct Bugs> data loaded Bugs> model compiled Bugs> model is initialized Bugs> 10000 updates took 0 s Bugs> monitor set Bugs> monitor set Bugs> monitor set Bugs> monitor set Bugs> 10000 updates took 0 s Bugs> mean sd val2.5pc median val97.5pc sample a 8428.0 890.5 6486.0 8418.0 10220.0 10000 b -407.4 74.07 -554.9 -407.6 -243.6 10000 Bugs> mean sd MC_error val2.5pc median val97.5pc start sample a 8428.0 890.5 56.85 6753.0 8425.0 10230.0 10001 10000 b -407.4 74.07 4.753 -554.6 -406.8 -266.0 10001 10000 Bugs> CODA files written Bugs> >The output is cars which is a matrix with 10 000 rows and 3 columns. The first column is the iteration number i.e. 10 001 to 20 000 (since the burn in period was 10 000, and the chain is run for a further 10 000 iterations). The second column is the values of the parameter a for each iteration, and the third column is the equivalent for the parameter b.
You can then proceed to manipulate the cars data set e.g.
apply(cars[,-1], 2, mean)gives the mean for a to be 8427.6630 and for b it is -407.3614 which agree with the direct output from LinBugs; and
plot(cars[,c(1,3)], type="l")which gives a plot of the trajectory of the MCMC for parameter b.