tests/checkmethods.R

In runjags: Interface Utilities, Model Templates, Parallel Computing Methods and Additional Distributions for MCMC Models in JAGS

# test all of the runjags dispatch methods with a toy example:

library(runjags)

runjags.options(inits.warning=FALSE, rng.warning=FALSE, blockignore.warning=FALSE, silent.jags=TRUE)

# Require for as.mcmc.list and niter:
library("coda")

model <- "model {
for(i in 1 : N){ #data# N
Y[i] ~ dnorm(true.y[i], precision); #data# Y
true.y[i] <- (m * X[i]) + c; #data# X
}
m ~ dunif(-1000,1000); #inits# m
c ~ dunif(-1000,1000);
precision ~ dexp(1);
#monitor# m, c, precision, dic, popt, pd
}"

model_parallel <- "model {
for(i in 1 : N){ #data# N
Y[i] ~ dnorm(true.y[i], precision); #data# Y
true.y[i] <- (m * X[i]) + c; #data# X
}
m ~ dunif(-1000,1000); #inits# m
c ~ dunif(-1000,1000);
precision ~ dexp(1);
#monitor# m, c, precision
}"

# Simulate the data
X <- 1:100
Y <- rnorm(length(X), 2*X + 10, 1)
N <- length(X)

initfunction <- function(chain) return(switch(chain, "1"=list(m=-10), "2"=list(m=10)))


# Get the JAGS path - rjags adds the JAGS path to the PATH in Windows...
try(library(rjags))
jagspath <- findjags()

# Only run the JAGS methods tests if we have found JAGS and have permission to run it:
if(jagspath!="JAGS not found" && testjags(jagspath)$JAGS.available){

	testnum <- 1
	cat('Testing the simple method\n')
	# Try the simple method and if it doesn't work give a warning but don't fail (likely to be permissions problems)
	success <- try({
		results <- run.jags(model, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='simple',temp=FALSE)
	})
	if(inherits(success, 'try-error')){
		cat("JAGS was found but the simple method failed; it is possible that there were permissions issues or similar.  Details as follows:\n")
		t <- testjags(silent=FALSE)
		cat(failed.jags('output')[[1]])
		print(file.info(jagspath))
		print(file.info(getwd())[,1:3])
		cat("All test methods (except possibly rjags) were skipped\n")
	}else{
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		stopifnot(niter(as.mcmc.list(results))==1000)

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results <- run.jags(model_parallel, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='parallel')
		stopifnot(niter(as.mcmc.list(results))==1000)

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results <- run.jags(model, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='interruptible')
		stopifnot(niter(as.mcmc.list(results))==1000)

		# Same as in checkinputs but it's for rjags there:
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results2 <- extend.jags(results, sample=1000, add.monitor="true.y", summarise=FALSE)
		stopifnot(nvar(as.mcmc.list(results2))==(3+N))
		stopifnot(varnames(as.mcmc.list(results2))[1]=='m' && varnames(as.mcmc.list(results2))[103]=='true.y[100]')

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results2 <- extend.jags(results, sample=1000, add.monitor=c("true.y", "dic"), summarise=FALSE)
		stopifnot(nvar(as.mcmc.list(results2))==(3+N))
		stopifnot(varnames(as.mcmc.list(results2))[1]=='m' && varnames(as.mcmc.list(results2))[103]=='true.y[100]')

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results2 <- extend.jags(results, sample=1000, add.monitor=c("true.y", "deviance"), summarise=FALSE)
		stopifnot(nvar(as.mcmc.list(results2))==(4+N))
		stopifnot(varnames(as.mcmc.list(results2))[1]=='m' && varnames(as.mcmc.list(results2))[104]=='deviance')

		# Same as in checkinputs but it's for rjags there:
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results2 <- extend.jags(results, sample=1000, add.monitor="true.y", summarise=FALSE, method='simple')
		stopifnot(nvar(as.mcmc.list(results2))==(3+N))
		stopifnot(varnames(as.mcmc.list(results2))[1]=='m' && varnames(as.mcmc.list(results2))[103]=='true.y[100]')

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results2 <- extend.jags(results, sample=1000, add.monitor=c("true.y", "dic"), summarise=FALSE)
		stopifnot(nvar(as.mcmc.list(results2))==(3+N))
		stopifnot(varnames(as.mcmc.list(results2))[1]=='m' && varnames(as.mcmc.list(results2))[103]=='true.y[100]')

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results2 <- extend.jags(results, sample=1000, add.monitor=c("true.y", "deviance"), summarise=FALSE)
		stopifnot(nvar(as.mcmc.list(results2))==(4+N))
		stopifnot(varnames(as.mcmc.list(results2))[1]=='m' && varnames(as.mcmc.list(results2))[104]=='deviance')

		# Snow gives problems here ... but it does work!
		#results <- run.jags(model, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='snow')
		#stopifnot(niter(as.mcmc.list(results))==1000)

		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		info <- run.jags(model, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='background')
		t <- 0
		repeat{
			# Change thin, chain and variables:
			s <- try(results <- results.jags(info, read.monitor='m', return.samples=100, recover.chains=2, summarise=FALSE))
			if(!inherits(s, 'try-error')) break
			if(t==30) stop("Timed out waiting for the bgparallel method")
			Sys.sleep(1)
			t <- t+1
		}
		stopifnot(niter(as.mcmc.list(results))==100)
		stopifnot(nvar(as.mcmc.list(results))==1)
		stopifnot(nchain(as.mcmc.list(results))==1)


		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		info <- run.jags(model_parallel, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='bgparallel')
		t <- 0
		repeat{
			s <- try(results <- results.jags(info))
			if(!inherits(s, 'try-error')) break
			if(t==30) stop("Timed out waiting for the bgparallel method")
			Sys.sleep(1)
			t <- t+1
		}
		stopifnot(niter(as.mcmc.list(results))==1000)

		# Check combine.mcmc does what it says on the tin:
		stopifnot(niter(combine.mcmc(results, return.samples=1000, collapse.chains=TRUE))==1000)
		stopifnot(niter(combine.mcmc(results, return.samples=11, collapse.chains=TRUE))==11)
		stopifnot(niter(combine.mcmc(results, return.samples=100, collapse.chains=FALSE))==100)
		stopifnot(niter(combine.mcmc(results, thin=1, collapse.chains=TRUE))==2000)
		stopifnot(niter(combine.mcmc(results, thin=10, collapse.chains=TRUE))==200)
		stopifnot((niter(combine.mcmc(results, thin=15, collapse.chains=TRUE))*15)>=2000)
		stopifnot(niter(combine.mcmc(results, thin=10, collapse.chains=FALSE))==100)

		# Check we can use the extend wrapper:
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		newres <- extend.jags(results, sample=0)
		stopifnot(newres$burnin==results$burnin)
		stopifnot(niter(as.mcmc.list(newres))==niter(as.mcmc.list(results)))

		# Check a single iteration works:
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results <- run.jags(model, n.chains=2, sample=1, burnin=1000, inits=initfunction, method='interruptible', summarise=FALSE)
		stopifnot(niter(as.mcmc.list(results))==1)

		# And that precision can be ignored:
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		results <- run.jags(model, monitor=c('m','c'), n.chains=2, sample=10, burnin=100, inits=initfunction, method='interruptible', summarise=FALSE, noread.monitor='precision')
		stopifnot(nvar(as.mcmc.list(results))==2)

		# Check the version number is correct:
		cat('Running method test number', testnum, '\n'); testnum <- testnum+1
		stopifnot(newres$runjags.version[1]==runjags:::runjagsprivate$runjagsversion)

	}
}else{
	cat("JAGS could not be called externally at the path: ", jagspath, "\n")
	cat("All test methods except possibly rjags and rjagsparallel were skipped\n")
}

testnum <- 1
if(require("rjags")){
	cat('Running rjags method test number', testnum, '\n'); testnum <- testnum+1
	results <- run.jags(model, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='rjags')
	stopifnot(niter(as.mcmc.list(results))==1000)

	cat('Running rjags method test number', testnum, '\n'); testnum <- testnum+1
	results <- run.jags(model, n.chains=2, sample=1, burnin=1000, inits=initfunction, method='rjags', summarise=FALSE)
	stopifnot(niter(as.mcmc.list(results))==1)

	cat('Running rjags method test number', testnum, '\n'); testnum <- testnum+1
	runjags.options(silent.jags=TRUE, silent.runjags=TRUE, debug=FALSE)
	output <- capture.output(results <- run.jags(model_parallel, n.chains=2, sample=1000, burnin=1000, inits=initfunction, method='rjparallel'))
	stopifnot(length(output)==0)
	stopifnot(niter(as.mcmc.list(results))==1000)
}else{
	cat("The rjags and rjagsparallel method checks were skipped as the rjags package is not installed\n")
}

cat("All methods checks passed\n")