From HPC

BUGS (Bayesian inference Using Gibbs Sampling) is a piece of computer software for the Bayesian analysis of complex statistical models using Markov chain Monte Carlo (MCMC) methods. It has been developing and maturing over the years, and is probably best known in its WinBugs incarnations. The latest version can run on Windows and Linux, as well as from inside the R statistical package.

Linbugs is Openbugs 3.0.3 classic interface.

You can run your Winbugs\Openbugs work on the HPC cluster using Linbugs which is the classic interface version and can be used from a terminal.

Openbugs web site http://mathstat.helsinki.fi/openbugs/

Openbugs manual http://mathstat.helsinki.fi/openbugs/Manuals/Manual.html

Contents 1 Opening Linbugs 2 Issues and Problems 3 Example 4 Using with R

Opening Linbugs

linbugs

Issues and Problems

The only issue know with the linux version is the command modelSaveLog('log.txt') cause the an error. This command should not be used, you will find that a file called run.out is create automatically which is the same as the log.

Example

Create the following files

data.txt

list(
  t = c( 9.43000E+01, 1.57000E+01, 6.29000E+01, 1.26000E+02, 
         5.24000E+00, 3.14000E+01, 1.05000E+00, 1.05000E+00, 
         2.10000E+00, 1.05000E+01 ), 
  x = c(  5, 1, 5, 14, 3, 19, 1, 1, 4, 22 ), 
  N = 10
)

model.txt

model
{
  for (i in 1 : N) {
    theta[i] ~ dgamma(alpha, beta)
    lambda[i] <- theta[i] * t[i]
    x[i] ~ dpois(lambda[i])
  }
  alpha ~ dexp(1)
  beta ~ dgamma(0.1, 1.0)
}

init1.txt

list(
  alpha = 1.00000E+00,
  beta = 1.00000E+00
)

Now open linbugs

linbugs

and enter the following commands (note, this script could be saved as a text file and run via the command 'linbugs < script.txt')

modelCheck('model.txt')
modelData('data.txt')
modelCompile(1)
modelInits('init1.txt', 1)
modelGenInits()
samplesBeg(501)
modelSetRN(314159)
modelUpdate(500, 1)
samplesSet(theta)
samplesSet(alpha)
samplesSet(beta)
samplesSet(deviance)
modelUpdate(500, 1)
samplesStats('*')
samplesCoda('*', 'coda')
modelQuit()

This will have create the following CODA files

codaCODAindex.txt
codaCODAchain1.txt

Using with R

This code uses two packages rbugs and coda. In theory rbugs should be enough, but it has a bug in reading in the coda files by it's self.

Create a r script file called 'pumps.r'

library(rbugs)
library(coda)

#Common variables

myWorkingDir = "/users/MyUserName/pumps"
myBugsWorkingDir = myWorkingDir
myN.chains = 1
myLinbugsBin = "/usr/local/openbugs/linbugs"

data(pumps)
pumps.data <- list(t = pumps$t, x = pumps$x, N = nrow(pumps))
pumps.model <- file.path(.path.package("rbugs"), "bugs/model", "pumps.bug")
pumps.inits <- file.path(.path.package("rbugs"), "bugs/inits", "pumps.txt")
inits <- list(dget(pumps.inits))
parameters <- c("theta", "alpha", "beta")
pumps.sim <- rbugs(data = pumps.data, 
 inits, 
 parameters,
 pumps.model, 
 n.chains = myN.chains, 
 n.iter = 1000,
 bugsWorkingDir = myBugsWorkingDir,
 workingDir = myWorkingDir,
 useWine = FALSE,
 linbugs = TRUE,
 bugs = myLinbugsBin,
 verbose = FALSE,
 genFilesOnly = TRUE) #genFilesOnly is important so we can run the bugs script manually and use coda library to read results
 
runBugs(bugs= myLinbugsBin, 
 script= paste(myBugsWorkingDir, "/script.txt", sep=""), 
 n.chains = myN.chains, 
 workingDir= myBugsWorkingDir, 
 linbugs = TRUE, 
 verbose = FALSE)

#Get results using the coda library
pumps.sim <- read.openbugs(stem= paste(myBugsWorkingDir, "/coda", sep=""), quiet=FALSE) 

str(pumps.sim)

This can then be run with the following command

R CMD BATCH pumps.r

This will create a file called pumps.r.Rout which should contain something similar to this at the end of the file

...
 ...Top part striped
...
> 
> pumps.sim <- read.openbugs(stem= paste(myBugsWorkingDir, "/coda", sep=""), quiet=FALSE) 
Abstracting alpha ... 499 valid values
Abstracting beta ... 499 valid values
Abstracting deviance ... 499 valid values
Abstracting theta[1] ... 499 valid values
Abstracting theta[2] ... 499 valid values
Abstracting theta[3] ... 499 valid values
Abstracting theta[4] ... 499 valid values
Abstracting theta[5] ... 499 valid values
Abstracting theta[6] ... 499 valid values
Abstracting theta[7] ... 499 valid values
Abstracting theta[8] ... 499 valid values
Abstracting theta[9] ... 499 valid values
Abstracting theta[10] ... 499 valid values
>  
> str(pumps.sim)
List of 1
 $ : mcmc [1:499, 1:13] 0.623 0.477 0.856 1.013 1.272 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:499] "502" "503" "504" "505" ...
  .. ..$ : chr [1:13] "alpha" "beta" "deviance" "theta[1]" ...
  ..- attr(*, "mcpar")= num [1:3] 502 1000 1
 - attr(*, "class")= chr "mcmc.list"
> 
> proc.time()
 user  system elapsed 
 0.532   0.096   0.600