Description Usage Arguments Details Value Note Author(s) References See Also Examples
View source: R/cudaMultireg.slice.R
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.
1 2 | cudaMultireg.slice(slicedata, ymaskdata, R, keep = 5, nu.e = 3,
fsave = NA, zprior=FALSE, rng = 0)
|
slicedata |
list(slice=slice, niislicets=niislicets, mask=mask, dsgn=dsgn);
input slice data used in simulation as returned by |
ymaskdata |
list(yn = yn, kin = kin, nreg = nreg);
masked and standardised slice data as returned by |
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: |
zprior |
boolean {T,F}; default {F} - use just a mean for |
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) |
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.
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 |
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.
Adelino Ferreira da Silva, Universidade Nova de Lisboa, Faculdade de Ciencias e Tecnologia, Portugal, afs@fct.unl.pt.
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.
read.fmrislice
,
read.Zsegslice
,
premask
,
pmeans.hcoef
,
regpostsim
,
plot.hcoef.post
,
post.simul.hist
,
post.ppm
,
post.tseries
,
post.randeff
,
post.shrinkage.mean
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.