req_suggested_packages <- c("rtdists", "microbenchmark", "reshape2", "ggplot2", "ggforce") pcheck <- lapply(req_suggested_packages, requireNamespace, quietly = TRUE) if (any(!unlist(pcheck))) { message("Required package(s) for this vignette are not available/installed and code will not be executed.") knitr::opts_chunk$set(eval = FALSE) }
knitr::opts_chunk$set( collapse = TRUE, error = TRUE, comment = "#>" )
Function pfddm()
evaluates the distribution function (or cumulative distribution function, CDF) for the Ratcliff diffusion decision model (DDM) using different methods for approximating the full CDF, which contains an infinite sum. This vignette provides an overview of the mathematical details of the different approximations implemented in pfddm()
as well as another approximation that is not included in pfddm()
for performance reasons. At the end of this vignette are timing benchmarks for the present approximation methods and comparison with existing approximation methods.
Our implementation of the DDM has the following parameters: $a \in (0, \infty)$ (threshold separation), $v \in (-\infty, \infty)$ (drift rate), $t_0 \in [0, \infty)$ (non-decision time/response time constant), $w \in (0, 1)$ (relative starting point), $sv \in (0, \infty)$ (inter-trial-variability of drift), and $\sigma \in (0, \infty)$ (diffusion coefficient of the underlying Wiener Process). Please note that for this vignette, we will refer to the inter-trial variability of drift as $\eta$ instead of $sv$ to make the notation in the equations less confusing.
There is one optional parameter in pfddm()
that can be used to indicate which method should be used in the function call: method
. For each approximation method we describe, we include the parameter settings for the function call to pfddm()
so that it uses the desired method. As this parameter is optional, leaving it blank results in the default method that is indicated later in this vignette. For general purpose use, we recommend ignoring this optional parameter so that the default settings are used as this will be the fastest and most stable algorithm for calculating the CDF.
Note that there are actually two cumulative distribution functions of the DDM: the CDF with respect to the upper boundary, and the CDF with respect to the lower boundary. Importantly, only the combined CDF of upper and lower boundary is a proper CDF in the sense that it reaches 1 at $t = \infty$. Following the precedent set by the literature, we use the general terminology "CDF" or "distribution function" to mean the defective cumulative distribution function with respect to the lower boundary (defective in the sense that it does not reach 1). Should the cumulative distribution function with respect to the upper boundary be required, it may be calculated using the simple transformation $F_\text{upper}(t ~|~ v, \eta, a, w) = F_\text{lower}(t ~|~ -v, \eta, a, 1-w)$.
Since the CDF of the DDM is widely used in parameter estimation usually involving numerical optimization (e.g., $\chi^2$ and Kolmogorov-Smirnov stastics), significant effort has been put into making its evaluation as fast as possible. However, the options for calculating the CDF are limited; one can either approximate an infinite sum or numerically solve a partial differential equation (PDE). The pfddm()
function uses the infinte sum method but not the PDE method; this vignette details these choices and demonstrates why the infinite sum method is preferable over the PDE method.
There are two main ways of calculating the CDF of the DDM: approximating the analytic form of CDF itself [@blurton2017first], and numerically solving a partial differential equation (PDE) whose solution is the CDF [@voss2008fast]. Like the probability density function (PDF) of the DDM, the analytic form of the CDF also contains an infinite sum (see the Math Vignette for more information on the PDF). We include the infinite sum method in pfddm()
but not the PDE method, and these reasons will be discussed in the following sections.
Both approximation methods inherently introduce error into the calculation (i.e., difference between the result of the approximated CDF and the "true" value of the CDF). Understanding this approximation error is important to minimizing the uncertainty in DDM parameter estimation (and other uses) that depend on calculating the CDF.
To use this method, set `method = "NCDF"` in the function call to `pfddm()`. We caution use of this approximation method, as the [default-method](#default-method) (using Mills ratio) is both faster and more stable.
To use this method, the user can ignore the optional parameter as this is the default method (internally, it sets `method = "Mills"` in the function call to `pfddm()`). This is the default method, and we recommend using it for the fastest and most stable algorithm for calculating the CDF of the DDM.
## PDE Method {#dist-pde} This section only provides a brief overview of the PDE method, highlighting how the PDE method may not be useful for our application; for more details on the PDE method, see @voss2008fast. Essentially, there exists a PDE whose solution is the CDF of the DDM; hence, solving the PDE allows one to calculate the CDF. The approximation error is controlled is via the step size in the numerical solution of the PDE. Compared to the infinite sum approach from [the previous section](#dist-inf), controlling the step size does not allow for directly bounding the approximation error. Whereas this method of controlling the approximation error may be suitable for some use cases, we suggest the infinite sum approach because it guarantees the precision of the result. In our testing with the implementation of the algorithm provided by @voss2008fast (using the function `rtdists::pdiffusion()` from the [`rtdists`](https://cran.r-project.org/package=rtdists) package), we found it to be inaccurate relative to the infinite sum method quite often. Increasing the `precision` parameter usually increased the accuracy of the PDE method (by reducing the step size in the numerical solver), but that resulted in very slow calculation times. It should be noted, however, that the method proposed by @voss2008fast is not only applicable to the case of inter-trial variability in the drift rate ($v$), but it can also be used in the cases of inter-trial variability in the starting point ($w$) and non-decision time ($t_0$).This method is unavailable in `fddm`. Instead, we recommend using the [default method](#default-method) (using Mills ratio) for calculating the CDF of the DDM by ignoring the optional parameter `method` in the function call to `pfddm()`. If the PDE method is desired, we recommend using the function `rtdists::pdiffusion()` from the [`rtdists`](https://cran.r-project.org/package=rtdists) package, as this function uses the PDE method to calculate the CDF.
We want to determine the performance of the two methods available in pfddm()
as well as comparing their performance to the implementations currently available in the literature: the R
Code provided by @blurton2017first, and the rtdists::pdiffusion()
function from the rtdists
package. Note that we only test implementations that include inter-trial variability in the drift rate (i.e., sv > 0
in pfddm()
).
Testing each implementation across a wide parameter space (defined in the code chunk below) will not only show the algorithms' overall viability for a practical set of parameters, but it will also allow for a more granular analysis of where each algorithm succeeds and struggles in the parameter space. To measure this viability, we will perform two benchmark tests: one where the response times are input as a vector to each function call, and one where the response times will be input individually to each function call.
# Define parameter space RT <- seq(0.1, 2, by = 0.1) A <- seq(0.5, 3.5, by = 0.5) V <- c(-5, -2, 0, 2, 5) t0 <- 0 W <- seq(0.3, 0.7, by = 0.1) SV <- c(0, 1, 2, 3.5) err_tol <- 1e-6 # this is the setting from rtdists
Note that each function call will assume that the response corresponds to the lower boundary of the DDM, and thus these benchmark tests only assess the performance of the "lower" distribution function. It is redundant to benchmark the "upper" distribution function because it is identical to the "lower" distribution function with the following mappings: $v \to -v$ and $w \to 1-w$; and these remapped values are already included in the parameter space (i.e., the values of $v$ are symmetric around $0$, and the values of $w$ are symmetric around $0.5$) .
For each combination of parameters in the parameter space, we run the microbenchmark::microbenchmark()
function from the package microbenchmark
1000 times for each implementation and only save the median benchmark time of these 1000. We will refer to these medians as simply the "benchmark times" for the remainder of this vignette. Running the benchmark tests this many times for each set of parameters generates a sufficient amount of benchmark data so that the median of these times is unlikely to be an outlier, either faster or slower than a typical run time.
First, we load the required packages.
# packages for generating benchmark data library("fddm") source(system.file("extdata", "Blurton_et_al_distribution.R", package = "fddm", mustWork = TRUE)) library("rtdists") library("microbenchmark") # packages for plotting library("reshape2") library("ggplot2") library("ggforce")
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.