R/RcppExports.R
In TasCL/ggdmc: Dynamic Model of Choice with Parallel Computation, and C++ Capabilities
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#' Map a parameter vector to an accumulator matrix
#'
#' An Rcpp function to get accumulator x parameter matrix for a design cell. For
#' two accumulator model, like drift-diffusion model, the accumulators are r1
#' and r2. For mulitple-accumulator model, like LBA model, the accumulators
#' are r1, r2, r3, etc. But because LBA probability density has yet
#' implemented, the function now only works for DDM. The part calling
#' \code{transform.dmc} is replaced by ordering the parameter sequence inside
#' \code{density.cpp}. The \code{n1order} switch is still kept, until the R's
#' \code{simulate.dmc} is upgraded to C++ version. This function is akin to
#' DMC's \code{p.df.dmc}. This is a Rcpp function.
#'
#' @param pVec a user supplied parameter vector or a sampler supplied theta/phi
#' vector.
#' @param cell a string indicating the name of a experimental condition
#' @param model an EAM model specification, created usually by model.dmc. Read
#' as a NumericVector, so attributes are kept.
#' @param n1order a boolean value indicating if swap accumulator row to use
#' "n1" order. This is for LBA model
#' @return an accumulator x EAM parameter matrix (without factor level).
#' @export
#' @examples
#' m1 <- model.dmc(
#' p.map = list(a="1",v="F",z="1",d="1",sz="1",sv="F",t0="1",st0="1"),
#' match.map = list(M=list(s1="r1", s2="r2")),
#' factors = list(S=c("s1","s2"), F=c("f1","f2")),
#' constants = c(st0=0,d=0),
#' responses = c("r1","r2"),
#' type = "rd")
#' # Parameter vector names are: ( see attr(,"p.vector") )
#' # [1] "a" "v.f1" "v.f2" "z" "sz" "sv.f1" "sv.f2" "t0"
#' #
#' # Constants are (see attr(,"constants") ):
#' # st0 d
#' # 0 0
#' #
#' # Model type = rd
#'
#' attr(m1,"p.vector")
#' # a v.f1 v.f2 z sz sv.f1 sv.f2 t0
#' # NA NA NA NA NA NA NA NA
#' pVec <- c(1,1,1,1, 1,1,1,1)
#' names(pVec) <- names(attr(m1, "p.vector"))
#'
#' accMat1 <- getAccumulatorMatrix(pVec, "s1.f1.r2", m1, FALSE)
#' ## col_names
#' ##row_names a v z d sz sv t0 st0
#' ## r1 1 1 0 0 1 1 1 0
#' ## r2 1 1 0 0 1 1 1 0
#' accMat2 <- getAccumulatorMatrix(pVec, "s1.f1.r2", m1, TRUE)
getAccumulatorMatrix <- function(pVec, cell, model, n1order = TRUE) {
.Call('ggdmc_getAccumulatorMatrix', PACKAGE = 'ggdmc', pVec, cell, model, n1order)
}
#' Compute Probability Density of Drift-Diffusion Model
#'
#' This function implements the equations in Voss, Rothermund, and Voss (2004).
#' These equations calculate Ratcliff's drift-diffusion model (RDM, 1978).
#' \code{ddmc} re-implements Voss, Rothermund, and Voss's (2004) equations A1
#' to A4 (page 1217) via Rcpp. This Rcpp function is akin to DMC's
#' \code{likelihood.dmc}. There is a \code{ddmc_parallel} with same arguments
#' using Open MPI to calculate RDM probability densities when data sets are
#' large (e.g., more than 5000 trials per condition) or when precision is set
#' high.
#'
#' @param x a data frame containing choices and RTs. This is akin to dnorm's x
#' argument.
#' @param pVec a parameter vector. For example,
#' p.vector <- c(a=1.25, v.f1=.20, v.f2=.33, z=.67, sz=.26, sv=.67, t0=.26)
#' @param precision a precision parameter. The larger, the higher precision to
#' estimate diffusion probability density. A general recommendation for the low
#' bound of precision is 2.5. If you need to use a precision higher than that,
#' you may want to consider using \code{ddmc_parallel}. A default value is set
#' at 2.5.
#' @param minLike a minimal log likelihood. If a estimated density is
#' less than minLike, the function returns minLike. A default value is set at
#' 1e-10.
#' @return a double vector with drift-diffusion probability density
#' @references Voss, A., Rothermund, K., & Voss, J. (2004). Interpreting the
#' parameters of the diffusion model: An empirical validation.
#' \emph{Memory & Cognition}, \bold{32(7)}, 1206-1220. \cr\cr
#' Ratcliff, R. (1978). A theory of memory retrival. \emph{Psychological
#' Review}, \bold{85}, 238-255.
#' @export
#' @examples
#' ## Set up a basic RDM design
#' m1 <- model.dmc(
#' p.map=list(a="1", v="1", z="1", d="1", sz="1", sv="1", t0="1", st0="1"),
#' constants=c(st0=0, d=0),
#' match.map=list(M=list(s1="r1", s2="r2")),
#' factors=list(S=c("s1", "s2")),
#' responses=c("r1", "r2"),
#' type="rd")
#'
#' pVec <- c(a=1, v=1, z=.5, sz=.25, sv=.2,t0=.15)
#'
#' ## Set up a model-data instance
#' raw.data <- simulate(m1, nsim=1e2, p.vector=pVec)
#' mdi <- data.model.dmc(raw.data, m1)
#'
#' ## Print probability densities on the screen
#' ddmc(mdi, pVec)
#' ## [1] 0.674039828 0.194299671 0.377715949 1.056043402 0.427508806
#' ## [6] 2.876521367 0.289301835 2.701315781 0.419909813 1.890380685
#' ## [11] 2.736803321 1.278188278 0.607926999 0.024858130 1.878312335
#' ## [16] 0.826923114 0.189952016 1.758620294 1.696875882 2.191290074
#' ## [21] 2.666750382 ...
#'
#' mdi.large <- data.model.dmc(simulate(m1, nsim=1e4, p.vector=pVec), m1)
#'
#' ## d.large <- dplyr::tbl_df(mdi.large)
#' ## d.large
#' ## Source: local data frame [20,000 x 3]
#' ##
#' ## S R RT
#' ## (fctr) (fctr) (dbl)
#' ## 1 s1 r2 0.2961630
#' ## 2 s1 r1 0.9745674
#' ## 3 s1 r2 0.6004945
#' ## 4 s1 r1 0.1937356
#' ## 5 s1 r1 0.3608844
#' ## 6 s1 r1 0.5343379
#' ## 7 s1 r1 0.2595707
#' ## 8 s1 r1 0.3247941
#' ## 9 s1 r1 0.4901325
#' ## 10 s1 r2 0.5011935
#' ## .. ... ... ...
#'
#' ## system.time(den1 <- ddmc(mdi.large, pVec))
#' ## user system elapsed
#' ## 0.028 0.004 0.046
#'
#' ## system.time(den2 <- ddmc_parallel(mdi.large, pVec))
#' ## user system elapsed
#' ## 0.200 0.000 0.023
#' ## all.equal(den1, den2)
#' ## [1] TRUE
ddmc <- function(x, pVec, precision = 2.5, minLike = 1e-10) {
.Call('ggdmc_ddmc', PACKAGE = 'ggdmc', x, pVec, precision, minLike)
}
#' @rdname ddmc
#' @export
ddmc_parallel <- function(x, pVec, precision = 2.5, minLike = 1e-10) {
.Call('ggdmc_ddmc_parallel', PACKAGE = 'ggdmc', x, pVec, precision, minLike)
}
#' Calculate Drift-diffusion Probability Density
#'
#' \code{g_minus} and \code{g_plus} implement A1 to A4 equations in Voss,
#' Rothermund, and Voss (2004). These equations calculate Ratcliff's
#' drift-diffusion model (1978). This source codes are derived from
#' Voss & Voss's fast-dm 30.2 in density.c.
#'
#' Two parallel functions \code{g_minus_parallel} and \code{g_plus_parallel},
#' using OpenMP libraries to do numerical integration. They resolve the
#' problem when high precision (> 10) is required.
#'
#' @param pVec a 9-element parameter (double) vector. The user has to follow
#' the sequence strictly. a, v, zr, d, sz, sv, t0, st0, RT, precision.
#' @references Voss, A., Rothermund, K., & Voss, J. (2004). Interpreting the
#' parameters of the diffusion model: A empirical validation
#' \emph{Memory and Cognition}, \bold{32(7)}, 1206--1220. \cr\cr
#' Ratcliff, R (1978). A theory of memory retrieval. \emph{Psychology Review},
#' \bold{85(2)}, 59--108.
#' @export
#' @examples
#' pVec1 <- c(a=2, v=2.5, zr=0.5, d=0, sz=0.3, sv=1, t0=0.3, st0=0,
#' RT=.550, precision=2.5)
#' g_minus(pVec1)
#' ## [1] 0.04965882
#' g_plus(pVec1)
#' ## [1] 2.146265
#' pVec2 <- c(a=2, v=2.5, zr=0.5, d=0, sz=0.3, sv=1, t0=0.3, st0=0,
#' RT=.550, precision=20)
#' ## system.time(d1 <- g_plus_parallel(pVec2))
#' ## user system elapsed
#' ## 135.708 0.000 12.022
#' ## system.time(d2 <- g_plus(pVec2))
#' ## user system elapsed
#' ## 104.516 0.000 104.490
g_minus <- function(pVec) {
.Call('ggdmc_g_minus', PACKAGE = 'ggdmc', pVec)
}
#' @rdname g_minus
#' @export
g_plus <- function(pVec) {
.Call('ggdmc_g_plus', PACKAGE = 'ggdmc', pVec)
}
#' @rdname g_minus
#' @export
g_minus_parallel <- function(pVec) {
.Call('ggdmc_g_minus_parallel', PACKAGE = 'ggdmc', pVec)
}
#' @rdname g_minus
#' @export
g_plus_parallel <- function(pVec) {
.Call('ggdmc_g_plus_parallel', PACKAGE = 'ggdmc', pVec)
}
#' Set up a DMC Sample for a Participant
#'
#' This is C++ function to initialse a DMC sample. The user usually should
#' not directly use it.
#'
#' @param nmc number of Markov Chain Monte Carlo iteration
#' @param pList prior distribution setting. This is R's p.prior
#' @param data a model data instance created by \code{data.model.dmc}
#' @param samples a DMC posterior sample/object
#' @param theta1 A user supplied initial theta cube
#' @param startPList A user supplied (different) prior distribution setting.
#' This is R's start.prior.
#' @param setting a list container to store all DMC setting
#' @export
initialise_data <- function(nmc, pList = NULL, data = NULL, samples = NULL, theta1 = NULL, startPList = NULL, setting = NULL) {
.Call('ggdmc_initialise_data', PACKAGE = 'ggdmc', nmc, pList, data, samples, theta1, startPList, setting)
}
#' Set up a DMC Sample for Multiple Participants
#'
#' This is C++ function to initialse a DMC sample with multiple participants.
#' The user usually should not directly use it.
#'
#' @param nmc number of Markov Chain Monte Carlo iteration
#' @param pList prior distribution setting. This is R's p.prior
#' @param data a model data instance created by \code{data.model.dmc}
#' @param samples a DMC posterior sample
#' @param theta1 A user supplied initial theta cube
#' @param startPList A user supplied (different) prior distribution setting.
#' This is R's start.prior.
#' @param phi1 A user supplied initial phi cube
#' @param ppList prior distribution setting, hyper level
#' @param hStartPList A user supplied (different) hyper-prior distribution
#' setting
#' @param setting a list container to store all DMC setting
#' @export
initialise_hyper <- function(nmc, pList = NULL, data = NULL, samples = NULL, theta1 = NULL, startPList = NULL, phi1 = NULL, ppList = NULL, hStartPList = NULL, setting = NULL) {
.Call('ggdmc_initialise_hyper', PACKAGE = 'ggdmc', nmc, pList, data, samples, theta1, startPList, phi1, ppList, hStartPList, setting)
}
#' Calculate Prior Probability Density for an EAM
#'
#' \code{dprior} is C++ function. It matches 5 string types: \code{tnorm},
#' \code{beta_lu}, \code{gamma_l}, \code{lnorm_l}, and \code{constant} to
#' determine which density functions to call (via R API). For truncated normal
#' density, \code{dprior} calls \code{dtn_scalar}, an internal Rcpp function
#' built specific for \pkg{ggdmc}. Whetehr log the probability density is
#' determined by the boolean \code{log} sent in via \code{p.prior}. This
#' function is akin to DMC's \code{log.prior.dmc}.
#'
#' @param pVec the user's supplied parameter vector or a sampler supplied
#' theta/phi vector.
#' @param pPrior a list of list usually created by prior.p.dmc to store the
#' prior parameter setting.
#' @return a double vector with probability density for each model parameter
#' @export
#' @examples
#' ## Use Drift-diffusion model as an example
#' ddm.prior <- prior.p.dmc(
#' dists = c("tnorm", "tnorm", "beta", "tnorm", "beta", "beta"),
#' p1 = c(a = 1, v = 0, z = 1, sz = 1, sv = 1, t0 = 1),
#' p2 = c(a = 1, v = 2, z = 1, sz = 1, sv = 1, t0 = 1),
#' lower = c(0,-5, NA, NA, 0, NA),
#' upper = c(2, 5, NA, NA, 2, NA))
#'
#' view(ddm.prior)
#' ## mean sd lower upper log dist untrans
#' ## a 1 1 0 2 1 tnorm identity
#' ## v 0 2 -5 5 1 tnorm identity
#' ## z 1 1 0 1 1 beta_lu identity
#' ## sz 1 1 -Inf Inf 1 tnorm identity
#' ## sv 1 1 0 2 1 beta_lu identity
#' ## t0 1 1 0 1 1 beta_lu identity
#'
#' ddm.pVec <- c(a=1.15, v=-0.10, z=0.74, sz=1.23, sv=0.11, t0=0.87)
#' dprior(ddm.pVec, ddm.prior)
#' ## a v z sz sv t0
#' ##-0.5484734 -1.6008386 0.0000000 -0.9453885 0.0000000 0.0000000
#'
#' ## Use LBA model as an example
#' lba.prior <- prior.p.dmc(
#' dists = c("tnorm", "tnorm", "tnorm", "tnorm", "tnorm", "tnorm"),
#' p1 = c(A=.4, B=.6, mean_v.true=1, mean_v.false=0, sd_v.true=.5, t0=.3),
#' p2 = c(A=.1, B=.1, mean_v.true=.2, mean_v.false=.2, sd_v.true=.1, t0=.05),
#' lower = c(0, 0, NA, NA, 0, .1),
#' upper = c(NA, NA, NA, NA, NA, 1))
#'
#' view(lba.prior)
#' ## mean sd lower upper log dist untrans
#' ## A 0.4 0.1 0 Inf 1 tnorm identity
#' ## B 0.6 0.1 0 Inf 1 tnorm identity
#' ## mean_v.true 1 0.2 -Inf Inf 1 tnorm identity
#' ## mean_v.false 0 0.2 -Inf Inf 1 tnorm identity
#' ## sd_v.true 0.5 0.1 0 Inf 1 tnorm identity
#' ## t0 0.3 0.05 0.1 1 1 tnorm identity
#'
#' lba.pVec <- c(A=0.398, B=0.614, mean_v.true=1.040,
#' mean_v.false=-0.032, sd_v.true=0.485, t0=0.271)
#' dprior(lba.pVec, lba.prior)
#' ## A B mean_v.true mean_v.false sd_v.true
#' ## 1.3834782 1.3738466 0.6704994 0.6776994 1.3723968
dprior <- function(pVec, pPrior) {
.Call('ggdmc_dprior', PACKAGE = 'ggdmc', pVec, pPrior)
}
#' Generate Random Numbers from Prior Probability Distribution
#'
#' \code{dprior} matches 5 types of string: \code{tnorm}, \code{beta_lu},
#' \code{gamma_l}, \code{lnorm_l}, and \code{constant} to determine which
#' density functions to call. \code{dprior} calls beta, gamma, log-normal
#' density functions via R API. For truncated normal density, \code{dprior}
#' calls \code{dtn_scalar}, an internal Rcpp function built specific for
#' ggdmc. Whetehr log the probability density is determined by the boolean
#' \code{log} sent in via p.prior. This function is akin to DMC's
#' \code{log.prior.dmc}.
#'
#' @param pPrior a p.prior list
#' @param n how many random number to generate
#' @return a double matrix with nrow equal to the number of random numbers and
#' ncol equal to the number of EAM parameters (npar)
#' @export
#' @examples
#' ddm.prior <- prior.p.dmc(
#' dists = c("tnorm", "tnorm", "beta", "tnorm", "beta", "beta"),
#' p1 = c(a = 1, v = 0, z = 1, sz = 1, sv = 1, t0 = 1),
#' p2 = c(a = 1, v = 2, z = 1, sz = 1, sv = 1, t0 = 1),
#' lower = c(0,-5, NA, NA, 0, NA),
#' upper = c(2, 5, NA, NA, 2, NA))
#'
#' view(ddm.prior)
#' ## mean sd lower upper log dist untrans
#' ## a 1 1 0 2 1 tnorm identity
#' ## v 0 2 -5 5 1 tnorm identity
#' ## z 1 1 0 1 1 beta_lu identity
#' ## sz 1 1 -Inf Inf 1 tnorm identity
#' ## sv 1 1 0 2 1 beta_lu identity
#' ## t0 1 1 0 1 1 beta_lu identity
#'
#' rprior(ddm.prior, 9)
#' ## a v z sz sv t0
#' ## [1,] 0.97413686 0.78446178 0.9975199 -0.5264946 0.5364492 0.55415052
#' ## [2,] 0.72870190 0.97151662 0.8516604 1.6008591 0.3399731 0.96520848
#' ## [3,] 1.63153685 1.96586939 0.9260939 0.7041254 0.4138329 0.78367440
#' ## [4,] 1.55866180 1.43657110 0.6152371 0.1290078 0.2957604 0.23027759
#' ## [5,] 1.32520281 -0.07328408 0.2051155 2.4040387 0.9663111 0.06127237
#' ## [6,] 0.49628528 -0.19374770 0.5142829 2.1452972 0.4335482 0.38410626
#' ## [7,] 0.03655549 0.77223432 0.1739831 1.4431507 0.6257398 0.63228368
#' ## [8,] 0.71197612 -1.15798082 0.8265523 0.3813370 0.4465184 0.23955415
#' ## [9,] 0.38049166 3.32132034 0.9888108 0.9684292 0.8437480 0.13502154
rprior <- function(pPrior, n) {
.Call('ggdmc_rprior', PACKAGE = 'ggdmc', pPrior, n)
}
#' Run a Bayesian EAM Model for Fixed-effect or Random-effect
#'
#' These are the C++ functions hierarchical Bayesian models. The user usually
#' should not call this function directly. Use \code{h.run.dmc} or
#' \code{run.dmc} wrapper instead. Two parallel versions of \code{run} uses
#' Open MPI to calculate diffusion probability density. Use them by setting
#' \code{cores} greater than 2 in \code{h.run.dmc}
#'
#' @param samples a DMC sample/object
#' @param setting a container to store all settings for running model fitting
#' @param debug a debugging switch to test chain randomisation. Default false
#' @return a DMC sample
#' @export
run_data <- function(samples, setting, debug = FALSE) {
.Call('ggdmc_run_data', PACKAGE = 'ggdmc', samples, setting, debug)
}
#' @rdname run_data
#' @export
run_data_parallel <- function(samples, setting, debug = FALSE) {
.Call('ggdmc_run_data_parallel', PACKAGE = 'ggdmc', samples, setting, debug)
}
#' @rdname run_data
#' @export
run_hyper <- function(samples, setting) {
.Call('ggdmc_run_hyper', PACKAGE = 'ggdmc', samples, setting)
}
#' @rdname run_data
#' @export
run_hyper_parallel <- function(samples, setting) {
.Call('ggdmc_run_hyper_parallel', PACKAGE = 'ggdmc', samples, setting)
}
assign_pp_pLists <- function(samples_s1, usePhi) {
.Call('ggdmc_assign_pp_pLists', PACKAGE = 'ggdmc', samples_s1, usePhi)
}
#' Sum and Log Probability Density of a EAM model
#'
#' Get log likelihood summed over a model data instance. The input data has
#' to be a data frame carrying a model specification, which is usually
#' created by \code{model.data.dmc} function.
#' \code{summed_log_likelihood_parallel} does calculation in parallel. Use it
#' when only the data set is big.
#'
#' @param pVec a parameter vector
#' @param data a model data instance
#' argument.
#' @return a double scalar
#' @export
#' @examples
#' m1 <- model.dmc(
#' p.map = list(a="1", v="1", z="1", d="1", sz="1", sv="1", t0="1", st0="1"),
#' constants = c(st0=0, d=0),
#' match.map = list(M=list(s1="r1", s2="r2")),
#' factors = list(S=c("s1", "s2")),
#' responses = c("r1", "r2"),
#' type = "rd")
#'
#' pVec <- c(a=1, v=1, z=0.5, sz=0.25, sv=0.2,t0=.15)
#'
#' ## Set up a model-data instance
#' dat <- simulate(m1, nsim=1e2, p.vector=pVec)
#' mdi <- data.model.dmc(dat, m1)
#' summed_log_likelihood(pVec, mdi)
#' ## [1] 0.3796048
summed_log_likelihood <- function(pVec, data) {
.Call('ggdmc_summed_log_likelihood', PACKAGE = 'ggdmc', pVec, data)
}
#' @rdname summed_log_likelihood
#' @export
summed_log_likelihood_parallel <- function(pVec, data) {
.Call('ggdmc_summed_log_likelihood_parallel', PACKAGE = 'ggdmc', pVec, data)
}
#' Sum and Log Prior Density of a EAM model
#'
#' Get log likelihood summed over all prior parameters
#'
#' @param pVec a parameter vector
#' @param pPrior p.prior
#' @export
#' @examples
#' pVec <- c(a=1, v=1, z=0.5, sz=0.25, sv=0.2,t0=.15)
#'
#' p.prior <- prior.p.dmc(
#' dists = rep("tnorm", 6),
#' p1 = c(a=2, v=2.5, z=0.5, sz=0.3, sv=1, t0=0.3),
#' p2 = c(a=0.5, v=.5, z=0.1, sz=0.1, sv=.3, t0=0.05) * 5,
#' lower = c(0,-5, 0, 0, 0, 0),
#' upper = c(5, 7, 2, 2, 2, 2))
#'
#' summed_log_prior(pVec, p.prior)
#' ## -3.224406
summed_log_prior <- function(pVec, pPrior) {
.Call('ggdmc_summed_log_prior', PACKAGE = 'ggdmc', pVec, pPrior)
}
TasCL/ggdmc documentation built on May 9, 2019, 4:19 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#' Map a parameter vector to an accumulator matrix
#'
#' An Rcpp function to get accumulator x parameter matrix for a design cell. For
#' two accumulator model, like drift-diffusion model, the accumulators are r1
#' and r2. For mulitple-accumulator model, like LBA model, the accumulators
#' are r1, r2, r3, etc. But because LBA probability density has yet
#' implemented, the function now only works for DDM. The part calling
#' \code{transform.dmc} is replaced by ordering the parameter sequence inside
#' \code{density.cpp}. The \code{n1order} switch is still kept, until the R's
#' \code{simulate.dmc} is upgraded to C++ version. This function is akin to
#' DMC's \code{p.df.dmc}. This is a Rcpp function.
#'
#' @param pVec a user supplied parameter vector or a sampler supplied theta/phi
#' vector.
#' @param cell a string indicating the name of a experimental condition
#' @param model an EAM model specification, created usually by model.dmc. Read
#' as a NumericVector, so attributes are kept.
#' @param n1order a boolean value indicating if swap accumulator row to use
#' "n1" order. This is for LBA model
#' @return an accumulator x EAM parameter matrix (without factor level).
#' @export
#' @examples
#' m1 <- model.dmc(
#' p.map = list(a="1",v="F",z="1",d="1",sz="1",sv="F",t0="1",st0="1"),
#' match.map = list(M=list(s1="r1", s2="r2")),
#' factors = list(S=c("s1","s2"), F=c("f1","f2")),
#' constants = c(st0=0,d=0),
#' responses = c("r1","r2"),
#' type = "rd")
#' # Parameter vector names are: ( see attr(,"p.vector") )
#' # [1] "a" "v.f1" "v.f2" "z" "sz" "sv.f1" "sv.f2" "t0"
#' #
#' # Constants are (see attr(,"constants") ):
#' # st0 d
#' # 0 0
#' #
#' # Model type = rd
#'
#' attr(m1,"p.vector")
#' # a v.f1 v.f2 z sz sv.f1 sv.f2 t0
#' # NA NA NA NA NA NA NA NA
#' pVec <- c(1,1,1,1, 1,1,1,1)
#' names(pVec) <- names(attr(m1, "p.vector"))
#'
#' accMat1 <- getAccumulatorMatrix(pVec, "s1.f1.r2", m1, FALSE)
#' ## col_names
#' ##row_names a v z d sz sv t0 st0
#' ## r1 1 1 0 0 1 1 1 0
#' ## r2 1 1 0 0 1 1 1 0
#' accMat2 <- getAccumulatorMatrix(pVec, "s1.f1.r2", m1, TRUE)
getAccumulatorMatrix <- function(pVec, cell, model, n1order = TRUE) {
.Call('ggdmc_getAccumulatorMatrix', PACKAGE = 'ggdmc', pVec, cell, model, n1order)
}
#' Compute Probability Density of Drift-Diffusion Model
#'
#' This function implements the equations in Voss, Rothermund, and Voss (2004).
#' These equations calculate Ratcliff's drift-diffusion model (RDM, 1978).
#' \code{ddmc} re-implements Voss, Rothermund, and Voss's (2004) equations A1
#' to A4 (page 1217) via Rcpp. This Rcpp function is akin to DMC's
#' \code{likelihood.dmc}. There is a \code{ddmc_parallel} with same arguments
#' using Open MPI to calculate RDM probability densities when data sets are
#' large (e.g., more than 5000 trials per condition) or when precision is set
#' high.
#'
#' @param x a data frame containing choices and RTs. This is akin to dnorm's x
#' argument.
#' @param pVec a parameter vector. For example,
#' p.vector <- c(a=1.25, v.f1=.20, v.f2=.33, z=.67, sz=.26, sv=.67, t0=.26)
#' @param precision a precision parameter. The larger, the higher precision to
#' estimate diffusion probability density. A general recommendation for the low
#' bound of precision is 2.5. If you need to use a precision higher than that,
#' you may want to consider using \code{ddmc_parallel}. A default value is set
#' at 2.5.
#' @param minLike a minimal log likelihood. If a estimated density is
#' less than minLike, the function returns minLike. A default value is set at
#' 1e-10.
#' @return a double vector with drift-diffusion probability density
#' @references Voss, A., Rothermund, K., & Voss, J. (2004). Interpreting the
#' parameters of the diffusion model: An empirical validation.
#' \emph{Memory & Cognition}, \bold{32(7)}, 1206-1220. \cr\cr
#' Ratcliff, R. (1978). A theory of memory retrival. \emph{Psychological
#' Review}, \bold{85}, 238-255.
#' @export
#' @examples
#' ## Set up a basic RDM design
#' m1 <- model.dmc(
#' p.map=list(a="1", v="1", z="1", d="1", sz="1", sv="1", t0="1", st0="1"),
#' constants=c(st0=0, d=0),
#' match.map=list(M=list(s1="r1", s2="r2")),
#' factors=list(S=c("s1", "s2")),
#' responses=c("r1", "r2"),
#' type="rd")
#'
#' pVec <- c(a=1, v=1, z=.5, sz=.25, sv=.2,t0=.15)
#'
#' ## Set up a model-data instance
#' raw.data <- simulate(m1, nsim=1e2, p.vector=pVec)
#' mdi <- data.model.dmc(raw.data, m1)
#'
#' ## Print probability densities on the screen
#' ddmc(mdi, pVec)
#' ## [1] 0.674039828 0.194299671 0.377715949 1.056043402 0.427508806
#' ## [6] 2.876521367 0.289301835 2.701315781 0.419909813 1.890380685
#' ## [11] 2.736803321 1.278188278 0.607926999 0.024858130 1.878312335
#' ## [16] 0.826923114 0.189952016 1.758620294 1.696875882 2.191290074
#' ## [21] 2.666750382 ...
#'
#' mdi.large <- data.model.dmc(simulate(m1, nsim=1e4, p.vector=pVec), m1)
#'
#' ## d.large <- dplyr::tbl_df(mdi.large)
#' ## d.large
#' ## Source: local data frame [20,000 x 3]
#' ##
#' ## S R RT
#' ## (fctr) (fctr) (dbl)
#' ## 1 s1 r2 0.2961630
#' ## 2 s1 r1 0.9745674
#' ## 3 s1 r2 0.6004945
#' ## 4 s1 r1 0.1937356
#' ## 5 s1 r1 0.3608844
#' ## 6 s1 r1 0.5343379
#' ## 7 s1 r1 0.2595707
#' ## 8 s1 r1 0.3247941
#' ## 9 s1 r1 0.4901325
#' ## 10 s1 r2 0.5011935
#' ## .. ... ... ...
#'
#' ## system.time(den1 <- ddmc(mdi.large, pVec))
#' ## user system elapsed
#' ## 0.028 0.004 0.046
#'
#' ## system.time(den2 <- ddmc_parallel(mdi.large, pVec))
#' ## user system elapsed
#' ## 0.200 0.000 0.023
#' ## all.equal(den1, den2)
#' ## [1] TRUE
ddmc <- function(x, pVec, precision = 2.5, minLike = 1e-10) {
.Call('ggdmc_ddmc', PACKAGE = 'ggdmc', x, pVec, precision, minLike)
}
#' @rdname ddmc
#' @export
ddmc_parallel <- function(x, pVec, precision = 2.5, minLike = 1e-10) {
.Call('ggdmc_ddmc_parallel', PACKAGE = 'ggdmc', x, pVec, precision, minLike)
}
#' Calculate Drift-diffusion Probability Density
#'
#' \code{g_minus} and \code{g_plus} implement A1 to A4 equations in Voss,
#' Rothermund, and Voss (2004). These equations calculate Ratcliff's
#' drift-diffusion model (1978). This source codes are derived from
#' Voss & Voss's fast-dm 30.2 in density.c.
#'
#' Two parallel functions \code{g_minus_parallel} and \code{g_plus_parallel},
#' using OpenMP libraries to do numerical integration. They resolve the
#' problem when high precision (> 10) is required.
#'
#' @param pVec a 9-element parameter (double) vector. The user has to follow
#' the sequence strictly. a, v, zr, d, sz, sv, t0, st0, RT, precision.
#' @references Voss, A., Rothermund, K., & Voss, J. (2004). Interpreting the
#' parameters of the diffusion model: A empirical validation
#' \emph{Memory and Cognition}, \bold{32(7)}, 1206--1220. \cr\cr
#' Ratcliff, R (1978). A theory of memory retrieval. \emph{Psychology Review},
#' \bold{85(2)}, 59--108.
#' @export
#' @examples
#' pVec1 <- c(a=2, v=2.5, zr=0.5, d=0, sz=0.3, sv=1, t0=0.3, st0=0,
#' RT=.550, precision=2.5)
#' g_minus(pVec1)
#' ## [1] 0.04965882
#' g_plus(pVec1)
#' ## [1] 2.146265
#' pVec2 <- c(a=2, v=2.5, zr=0.5, d=0, sz=0.3, sv=1, t0=0.3, st0=0,
#' RT=.550, precision=20)
#' ## system.time(d1 <- g_plus_parallel(pVec2))
#' ## user system elapsed
#' ## 135.708 0.000 12.022
#' ## system.time(d2 <- g_plus(pVec2))
#' ## user system elapsed
#' ## 104.516 0.000 104.490
g_minus <- function(pVec) {
.Call('ggdmc_g_minus', PACKAGE = 'ggdmc', pVec)
}
#' @rdname g_minus
#' @export
g_plus <- function(pVec) {
.Call('ggdmc_g_plus', PACKAGE = 'ggdmc', pVec)
}
#' @rdname g_minus
#' @export
g_minus_parallel <- function(pVec) {
.Call('ggdmc_g_minus_parallel', PACKAGE = 'ggdmc', pVec)
}
#' @rdname g_minus
#' @export
g_plus_parallel <- function(pVec) {
.Call('ggdmc_g_plus_parallel', PACKAGE = 'ggdmc', pVec)
}
#' Set up a DMC Sample for a Participant
#'
#' This is C++ function to initialse a DMC sample. The user usually should
#' not directly use it.
#'
#' @param nmc number of Markov Chain Monte Carlo iteration
#' @param pList prior distribution setting. This is R's p.prior
#' @param data a model data instance created by \code{data.model.dmc}
#' @param samples a DMC posterior sample/object
#' @param theta1 A user supplied initial theta cube
#' @param startPList A user supplied (different) prior distribution setting.
#' This is R's start.prior.
#' @param setting a list container to store all DMC setting
#' @export
initialise_data <- function(nmc, pList = NULL, data = NULL, samples = NULL, theta1 = NULL, startPList = NULL, setting = NULL) {
.Call('ggdmc_initialise_data', PACKAGE = 'ggdmc', nmc, pList, data, samples, theta1, startPList, setting)
}
#' Set up a DMC Sample for Multiple Participants
#'
#' This is C++ function to initialse a DMC sample with multiple participants.
#' The user usually should not directly use it.
#'
#' @param nmc number of Markov Chain Monte Carlo iteration
#' @param pList prior distribution setting. This is R's p.prior
#' @param data a model data instance created by \code{data.model.dmc}
#' @param samples a DMC posterior sample
#' @param theta1 A user supplied initial theta cube
#' @param startPList A user supplied (different) prior distribution setting.
#' This is R's start.prior.
#' @param phi1 A user supplied initial phi cube
#' @param ppList prior distribution setting, hyper level
#' @param hStartPList A user supplied (different) hyper-prior distribution
#' setting
#' @param setting a list container to store all DMC setting
#' @export
initialise_hyper <- function(nmc, pList = NULL, data = NULL, samples = NULL, theta1 = NULL, startPList = NULL, phi1 = NULL, ppList = NULL, hStartPList = NULL, setting = NULL) {
.Call('ggdmc_initialise_hyper', PACKAGE = 'ggdmc', nmc, pList, data, samples, theta1, startPList, phi1, ppList, hStartPList, setting)
}
#' Calculate Prior Probability Density for an EAM
#'
#' \code{dprior} is C++ function. It matches 5 string types: \code{tnorm},
#' \code{beta_lu}, \code{gamma_l}, \code{lnorm_l}, and \code{constant} to
#' determine which density functions to call (via R API). For truncated normal
#' density, \code{dprior} calls \code{dtn_scalar}, an internal Rcpp function
#' built specific for \pkg{ggdmc}. Whetehr log the probability density is
#' determined by the boolean \code{log} sent in via \code{p.prior}. This
#' function is akin to DMC's \code{log.prior.dmc}.
#'
#' @param pVec the user's supplied parameter vector or a sampler supplied
#' theta/phi vector.
#' @param pPrior a list of list usually created by prior.p.dmc to store the
#' prior parameter setting.
#' @return a double vector with probability density for each model parameter
#' @export
#' @examples
#' ## Use Drift-diffusion model as an example
#' ddm.prior <- prior.p.dmc(
#' dists = c("tnorm", "tnorm", "beta", "tnorm", "beta", "beta"),
#' p1 = c(a = 1, v = 0, z = 1, sz = 1, sv = 1, t0 = 1),
#' p2 = c(a = 1, v = 2, z = 1, sz = 1, sv = 1, t0 = 1),
#' lower = c(0,-5, NA, NA, 0, NA),
#' upper = c(2, 5, NA, NA, 2, NA))
#'
#' view(ddm.prior)
#' ## mean sd lower upper log dist untrans
#' ## a 1 1 0 2 1 tnorm identity
#' ## v 0 2 -5 5 1 tnorm identity
#' ## z 1 1 0 1 1 beta_lu identity
#' ## sz 1 1 -Inf Inf 1 tnorm identity
#' ## sv 1 1 0 2 1 beta_lu identity
#' ## t0 1 1 0 1 1 beta_lu identity
#'
#' ddm.pVec <- c(a=1.15, v=-0.10, z=0.74, sz=1.23, sv=0.11, t0=0.87)
#' dprior(ddm.pVec, ddm.prior)
#' ## a v z sz sv t0
#' ##-0.5484734 -1.6008386 0.0000000 -0.9453885 0.0000000 0.0000000
#'
#' ## Use LBA model as an example
#' lba.prior <- prior.p.dmc(
#' dists = c("tnorm", "tnorm", "tnorm", "tnorm", "tnorm", "tnorm"),
#' p1 = c(A=.4, B=.6, mean_v.true=1, mean_v.false=0, sd_v.true=.5, t0=.3),
#' p2 = c(A=.1, B=.1, mean_v.true=.2, mean_v.false=.2, sd_v.true=.1, t0=.05),
#' lower = c(0, 0, NA, NA, 0, .1),
#' upper = c(NA, NA, NA, NA, NA, 1))
#'
#' view(lba.prior)
#' ## mean sd lower upper log dist untrans
#' ## A 0.4 0.1 0 Inf 1 tnorm identity
#' ## B 0.6 0.1 0 Inf 1 tnorm identity
#' ## mean_v.true 1 0.2 -Inf Inf 1 tnorm identity
#' ## mean_v.false 0 0.2 -Inf Inf 1 tnorm identity
#' ## sd_v.true 0.5 0.1 0 Inf 1 tnorm identity
#' ## t0 0.3 0.05 0.1 1 1 tnorm identity
#'
#' lba.pVec <- c(A=0.398, B=0.614, mean_v.true=1.040,
#' mean_v.false=-0.032, sd_v.true=0.485, t0=0.271)
#' dprior(lba.pVec, lba.prior)
#' ## A B mean_v.true mean_v.false sd_v.true
#' ## 1.3834782 1.3738466 0.6704994 0.6776994 1.3723968
dprior <- function(pVec, pPrior) {
.Call('ggdmc_dprior', PACKAGE = 'ggdmc', pVec, pPrior)
}
#' Generate Random Numbers from Prior Probability Distribution
#'
#' \code{dprior} matches 5 types of string: \code{tnorm}, \code{beta_lu},
#' \code{gamma_l}, \code{lnorm_l}, and \code{constant} to determine which
#' density functions to call. \code{dprior} calls beta, gamma, log-normal
#' density functions via R API. For truncated normal density, \code{dprior}
#' calls \code{dtn_scalar}, an internal Rcpp function built specific for
#' ggdmc. Whetehr log the probability density is determined by the boolean
#' \code{log} sent in via p.prior. This function is akin to DMC's
#' \code{log.prior.dmc}.
#'
#' @param pPrior a p.prior list
#' @param n how many random number to generate
#' @return a double matrix with nrow equal to the number of random numbers and
#' ncol equal to the number of EAM parameters (npar)
#' @export
#' @examples
#' ddm.prior <- prior.p.dmc(
#' dists = c("tnorm", "tnorm", "beta", "tnorm", "beta", "beta"),
#' p1 = c(a = 1, v = 0, z = 1, sz = 1, sv = 1, t0 = 1),
#' p2 = c(a = 1, v = 2, z = 1, sz = 1, sv = 1, t0 = 1),
#' lower = c(0,-5, NA, NA, 0, NA),
#' upper = c(2, 5, NA, NA, 2, NA))
#'
#' view(ddm.prior)
#' ## mean sd lower upper log dist untrans
#' ## a 1 1 0 2 1 tnorm identity
#' ## v 0 2 -5 5 1 tnorm identity
#' ## z 1 1 0 1 1 beta_lu identity
#' ## sz 1 1 -Inf Inf 1 tnorm identity
#' ## sv 1 1 0 2 1 beta_lu identity
#' ## t0 1 1 0 1 1 beta_lu identity
#'
#' rprior(ddm.prior, 9)
#' ## a v z sz sv t0
#' ## [1,] 0.97413686 0.78446178 0.9975199 -0.5264946 0.5364492 0.55415052
#' ## [2,] 0.72870190 0.97151662 0.8516604 1.6008591 0.3399731 0.96520848
#' ## [3,] 1.63153685 1.96586939 0.9260939 0.7041254 0.4138329 0.78367440
#' ## [4,] 1.55866180 1.43657110 0.6152371 0.1290078 0.2957604 0.23027759
#' ## [5,] 1.32520281 -0.07328408 0.2051155 2.4040387 0.9663111 0.06127237
#' ## [6,] 0.49628528 -0.19374770 0.5142829 2.1452972 0.4335482 0.38410626
#' ## [7,] 0.03655549 0.77223432 0.1739831 1.4431507 0.6257398 0.63228368
#' ## [8,] 0.71197612 -1.15798082 0.8265523 0.3813370 0.4465184 0.23955415
#' ## [9,] 0.38049166 3.32132034 0.9888108 0.9684292 0.8437480 0.13502154
rprior <- function(pPrior, n) {
.Call('ggdmc_rprior', PACKAGE = 'ggdmc', pPrior, n)
}
#' Run a Bayesian EAM Model for Fixed-effect or Random-effect
#'
#' These are the C++ functions hierarchical Bayesian models. The user usually
#' should not call this function directly. Use \code{h.run.dmc} or
#' \code{run.dmc} wrapper instead. Two parallel versions of \code{run} uses
#' Open MPI to calculate diffusion probability density. Use them by setting
#' \code{cores} greater than 2 in \code{h.run.dmc}
#'
#' @param samples a DMC sample/object
#' @param setting a container to store all settings for running model fitting
#' @param debug a debugging switch to test chain randomisation. Default false
#' @return a DMC sample
#' @export
run_data <- function(samples, setting, debug = FALSE) {
.Call('ggdmc_run_data', PACKAGE = 'ggdmc', samples, setting, debug)
}
#' @rdname run_data
#' @export
run_data_parallel <- function(samples, setting, debug = FALSE) {
.Call('ggdmc_run_data_parallel', PACKAGE = 'ggdmc', samples, setting, debug)
}
#' @rdname run_data
#' @export
run_hyper <- function(samples, setting) {
.Call('ggdmc_run_hyper', PACKAGE = 'ggdmc', samples, setting)
}
#' @rdname run_data
#' @export
run_hyper_parallel <- function(samples, setting) {
.Call('ggdmc_run_hyper_parallel', PACKAGE = 'ggdmc', samples, setting)
}
assign_pp_pLists <- function(samples_s1, usePhi) {
.Call('ggdmc_assign_pp_pLists', PACKAGE = 'ggdmc', samples_s1, usePhi)
}
#' Sum and Log Probability Density of a EAM model
#'
#' Get log likelihood summed over a model data instance. The input data has
#' to be a data frame carrying a model specification, which is usually
#' created by \code{model.data.dmc} function.
#' \code{summed_log_likelihood_parallel} does calculation in parallel. Use it
#' when only the data set is big.
#'
#' @param pVec a parameter vector
#' @param data a model data instance
#' argument.
#' @return a double scalar
#' @export
#' @examples
#' m1 <- model.dmc(
#' p.map = list(a="1", v="1", z="1", d="1", sz="1", sv="1", t0="1", st0="1"),
#' constants = c(st0=0, d=0),
#' match.map = list(M=list(s1="r1", s2="r2")),
#' factors = list(S=c("s1", "s2")),
#' responses = c("r1", "r2"),
#' type = "rd")
#'
#' pVec <- c(a=1, v=1, z=0.5, sz=0.25, sv=0.2,t0=.15)
#'
#' ## Set up a model-data instance
#' dat <- simulate(m1, nsim=1e2, p.vector=pVec)
#' mdi <- data.model.dmc(dat, m1)
#' summed_log_likelihood(pVec, mdi)
#' ## [1] 0.3796048
summed_log_likelihood <- function(pVec, data) {
.Call('ggdmc_summed_log_likelihood', PACKAGE = 'ggdmc', pVec, data)
}
#' @rdname summed_log_likelihood
#' @export
summed_log_likelihood_parallel <- function(pVec, data) {
.Call('ggdmc_summed_log_likelihood_parallel', PACKAGE = 'ggdmc', pVec, data)
}
#' Sum and Log Prior Density of a EAM model
#'
#' Get log likelihood summed over all prior parameters
#'
#' @param pVec a parameter vector
#' @param pPrior p.prior
#' @export
#' @examples
#' pVec <- c(a=1, v=1, z=0.5, sz=0.25, sv=0.2,t0=.15)
#'
#' p.prior <- prior.p.dmc(
#' dists = rep("tnorm", 6),
#' p1 = c(a=2, v=2.5, z=0.5, sz=0.3, sv=1, t0=0.3),
#' p2 = c(a=0.5, v=.5, z=0.1, sz=0.1, sv=.3, t0=0.05) * 5,
#' lower = c(0,-5, 0, 0, 0, 0),
#' upper = c(5, 7, 2, 2, 2, 2))
#'
#' summed_log_prior(pVec, p.prior)
#' ## -3.224406
summed_log_prior <- function(pVec, pPrior) {
.Call('ggdmc_summed_log_prior', PACKAGE = 'ggdmc', pVec, pPrior)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.