cudaMultireg.slice: CUDA Parallel Implementation of a Bayesian Multilevel Model...

In cudaBayesreg: CUDA Parallel Implementation of a Bayesian Multilevel Model for fMRI Data Analysis

Description

cudaMultireg.slice provides an interface to a CUDA implementation of a Bayesian multilevel model for the analysis of brain fMRI data. cudaMultireg.slice processes a single slice.

Usage

cudaMultireg.slice(slicedata, ymaskdata, R, keep = 5, nu.e = 3,
	fsave = NA, zprior=FALSE, rng = 0)

Arguments

slicedata	list(slice=slice, niislicets=niislicets, mask=mask, dsgn=dsgn); input slice data used in simulation as returned by read.fmrislice. See read.fmrislice for indications on how to process user defined datasets.
ymaskdata	list(yn = yn, kin = kin, nreg = nreg); masked and standardised slice data as returned by premask
R	number of MCMC draws
keep	MCMC thinning parameter: keep every keepth draw (def: 5)
nu.e	d.f. parameter for regression error variance prior (def: 3)
fsave	filename for saving the MCMC simulation (def: NULL do not save)
zprior	boolean {T,F}; default {F} - use just a mean for Z
rng	integer {0,1,2}: type of RNG to use {0-Marsaglia Multicarry, 1-R. P. Brent xorgens, 2-Mersenne Twister MT19937-64}; (def. 0-Marsaglia Multicarry)

Details

The statistical model implemented in CUDA was specified as a Gibbs Sampler for hierarchical linear models with a normal prior. This model has been analysed by Rossi, Allenby and McCulloch in Bayesian Statistics and Marketing, Chapt. 3, and is referred to as rhierLinearModel in the R package bayesm. The main computational work is done in parallel on a CUDA capable GPU. Each thread is responsible for fitting a general linear model at each voxel. The CUDA implementation has the following system requirements: nvcc NVIDIA Cuda Compiler driver, g++ GNU compiler (nvcc compatible version). The package includes source code files to build the library ‘newmat11.so’. This is a matrix library by R. B. Davies used in the package's host C++ code. The package includes three optional CUDA-based RNGs. Marsaglia's multicarry RNG follows the R implementation, is the fastest one, and is used by default; Brent's RNG has higher quality but is not-so-fast; Matsumoto's Mersenne Twister is slow. The data sets used in the examples are available in the R package cudaBayesregData.

Value

a list containing

betadraw	nreg x nvar x R/keep array of individual regression coef draws
taudraw	R/keep x nreg array of error variance draws
Deltadraw	R/keep x nz x nvar array of Deltadraws
Vbetadraw	R/keep x nvar*nvar array of Vbeta draws

Note

The statistical model may be specified as follows.

Model: length(regdata) regression equations.
y_i = X_ibeta_i + e_i. e_i ~ N(0,tau_i). nvar X vars in each equation.

Priors:
tau_i ~ nu.e*ssq_i/χ^2_{nu.e}. tau_i is the variance of e_i.
beta_i ~ N(ZDelta[i,],V_{beta}).
Note: ZDelta is the matrix Z * Delta; [i,] refers to ith row of this product.

vec(Delta) given V_{beta} ~ N(vec(Deltabar),V_{beta} (x) A^{-1}).
V_{beta} ~ IW(nu,V).
Delta, Deltabar are nz x nvar. A is nz x nz. V_{beta} is nvar x nvar.

By default we suppose that we don't have any z vars, Z=iota (nreg x 1).

Simulated objects are specified as in bayesm with classes bayesm.mat and bayesm.var. S3 methods to summarize marginal distributions given an array of draws are then compatible with those of bayesm (see Examples).
Summaries will be invoked by a call to the generic summary function as in summary(object) where object is of class bayesm.mat or bayesm.var.
A new S3 method (hcoef.post) is specified for dispatching betadraw plots.

Author(s)

Adelino Ferreira da Silva, Universidade Nova de Lisboa, Faculdade de Ciencias e Tecnologia, Portugal, [e-mail].

References

Adelino R. Ferreira da Silva (2011). “cudaBayesreg: Parallel Implementation of a Bayesian Multilevel Model for fMRI Data Analysis.” Journal of Statistical Software, 44(4), 1–24. URL http://www.jstatsoft.org/v44/i04/.

Adelino Ferreira da Silva (2011). cudaBayesregData: Data sets for the examples used in the package cudaBayesreg, R package version 0.3-10. URL http://CRAN.R-project.org/package=cudaBayesregData.

Adelino Ferreira da Silva (2011). “A Bayesian Multilevel Model for fMRI Data Analysis.”, Computer Methods and Programs in Biomedicine, 102,(3), 238–252.

Adelino Ferreira da Silva (2010). “cudaBayesreg: Bayesian Computation in CUDA.”, The R Journal, 2/2, 48-55. URL http://journal.r-project.org/archive/2010-2/RJournal_2010-2_Ferreira~da~Silva.pdf.

Rossi, Allenby and McCulloch. Bayesian Statistics and Marketing, Chapter 3. URL http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html.

Davies, R.B. (1994). Writing a matrix package in C++. In OON-SKI'94: The second annual object-oriented numerics conference, pp 207-213. Rogue Wave Software, Corvallis. URL http://www.robertnz.net/cpp\_site.html.

Richard. P. Brent. Some long-period random number generators using shifts and xors, Preprint: 2 July 2007.

Brandon Whitcher, Volker Schmid and Andrew Thornton (2011). oro.nifti: Rigorous - NIfTI Input / Output, R package version 0.2.5. URL http://CRAN.R-project.org/package=oro.nifti.

See Also

read.fmrislice, read.Zsegslice, premask, pmeans.hcoef, regpostsim, plot.hcoef.post, post.simul.hist, post.ppm, post.tseries, post.randeff, post.shrinkage.mean

Examples

## Not run: 
## Simulation using the visual/auditory test dataset "fmri"  
slicedata <- read.fmrislice(fbase="fmri", slice=3, swap=FALSE)
ymaskdata <- premask(slicedata)
fsave <- paste(tempdir(),"/simultest1",fileext = ".sav", sep="")
out <- cudaMultireg.slice(slicedata, ymaskdata, R=2000, keep=5, nu.e=3,
  fsave=fsave, zprior=FALSE, rng=0 )
## Post-processing simulation
post.ppm(out=out, slicedata=slicedata, ymaskdata=ymaskdata, vreg=2)
post.ppm(out=out, slicedata=slicedata, ymaskdata=ymaskdata, vreg=4)
## "bayesm" summaries 
require("bayesm")
summary(out$betadraw)
summary(out$Deltadraw)
plot(out$Deltadraw)
summary(out$Vbetadraw)
##
## Random effects simulation using the SPM auditory dataset "swrfM*" 
fbase <- "swrfM"
slice <- 21
slicedata <- read.fmrislice(fbase=fbase, slice=slice, swap=FALSE )
ymaskdata <- premask(slicedata)
fsave <- paste(tempdir(),"/simultest3",fileext = ".sav", sep="")
out <- cudaMultireg.slice(slicedata, ymaskdata, R=2000, keep=5, nu.e=3,
  fsave=fsave, zprior=TRUE, rng=1)
post.ppm(out=out, slicedata=slicedata, ymaskdata=ymaskdata, vreg=2)

## End(Not run)

cudaBayesreg documentation built on May 29, 2017, 6:19 p.m.