cudaBayesreg (version 0.3-6)

cudaMultireg.slice: 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.

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

ymaskdata

list(yn = yn, kin = kin, nreg = nreg); masked and standardized 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)

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

concept

• bayes
• MCMC
• Gibbs Sampling
• hierarchical models
• linear model

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.

References

Adelino Ferreira da Silva, A Bayesian Multilevel Model for fMRI Data Analysis, Computer Methods and Programs in Biomedicine, to be published. Rossi, Allenby and McCulloch, Bayesian Statistics and Marketing, Chapter 3. 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. http://www.robertnz.net/cpp_site.html. Richard. P. Brent, Some long-period random number generators using shifts and xors, Preprint: 2 July 2007.

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

## Simulation using the visual/auditory test dataset "fmri"  
slicedata <- read.fmrislice(fbase="fmri", slice=3, swap=FALSE)
ymaskdata <- premask(slicedata)
fsave <- "/tmp/simultest1.sav"
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 <- "/tmp/simultest3.sav"
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)