Function dfddm evaluates the density function (or probability density function, PDF) for the Ratcliff diffusion decision model (DDM) using different methods for approximating the full PDF, which contains an infinite sum. An empirical validation of the implemented methods is provided in the Validity Vignette. Timing benchmarks for the present methods and comparison with existing methods are provided in the Benchmark Vignette. Two examples of using dfddm for parameter estimation are provided in the Example Vignette.
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), and \(sv \in (0, \infty)\) (inter-trial-variability of drift).
There are several different methods for approximating the PDF of the DDM, and there are three optional parameters in dfddm that can be used to indicate which method should be in the function call: scale, n_terms_small, and summation_small. For each method we describe, we include the parameter settings for the function call to dfddm so that it uses the desired method. As these parameters are optional, leaving them blank results in the default method that is indicated later in this vignette. For general purpose use, we recommend ignoring these optional parameters so that the default settings are used as this will be the fastest algorithm. Note that precedence for the optional parameters is first given to checking if the default method is selected. If not, precedence is then given to the scale parameter value; for example, scale = "large" will ignore the summation_small input value.
Since the DDM is widely used in parameter estimation usually involving numerical optimization, significant effort has been put into making the evaluation of its density as fast as possible. However, the density function for the DDM is notorious for containing an unavoidable infinite sum; hence, the literature has produced a few different methods of approximating the density. This vignette details the various methods used in the literature to approximate the infinite sum in the density function for the DDM.
In the seminal book where the density function originates, Feller (1968) explains the derivation from first principles. In this derivation there is a step that requires taking a limit, and Feller (1968) provides two different – but equivalent – limiting processes that yield two different – but equal – forms of the density function. Both of these forms contain an infinite sum, and they are known as the large-time approximation and the small-time approximation because the former is on average faster when calculating the density for large response times and the latter is on average faster when calculating the density for small response times.
When the drift rate is held constant (i.e. \(sv = 0\)), the density function for the DDM is often written in a factorized form (Navarro and Fuss 2009): \[\begin{equation} f(t ~|~ v, a, w) = \frac{1}{a^2} \exp \left( -vaw -\frac{v^2 t}{2} \right) f_i \left( \frac{t}{a^2} ~\Big\vert~ 0, 1, w \right), \nonumber \end{equation}\] where \(f_i(\frac{t}{a^2} | 0, 1, w)\) determines whether the large-time or small-time model will be used: \[\begin{equation} \begin{aligned} f_{i=\ell} (t ~|~ 0, 1, w) &= \pi \sum_{j = 1}^{\infty} j \exp \left( -\frac{j^2 \pi^2 t}{2a^2} \right) \sin \left( j w \pi \right),\\ f_{i=s} (t ~|~ 0, 1, w) &= \frac{1}{\sqrt{2 \pi t^3}} \sum_{j = -\infty}^{\infty} (w + 2j) \exp \left( -\frac{(w + 2j)^2}{2t} \right). \end{aligned} \nonumber \end{equation}\]
In an effort to simplify the terms inside the infinite summations as much as possible, we instead rewrite the constant drift rate density function as two separate functions without the factorization: \[\begin{align} f_\ell(t | v, a, w) &= \frac{\pi}{a^2} e^{ \left( -vaw-\frac{v^2 t}{2} \right)} \sum_{j = 1}^{\infty} j \sin \left( j w \pi \right) \exp{ \left( -\frac{j^2 \pi^2 t}{2a^2} \right)}, \label{eq:con-l} \\ f_s(t | v, a, w) &= \frac{a}{\sqrt{2 \pi t^3}} e^{ \left( -vaw-\frac{v^2 t}{2} \right)} \sum_{j = -\infty}^{\infty} (w + 2j) \exp{ \left( -\frac{a^2}{2t} \left( w + 2j \right)^2 \right)}. \label{eq:con-s} \end{align}\]
In addition to having large-time and small-time variants, there exist two mathematically equivalent formulations for the infinite summation in the small-time density functions. The details and proof of equivalence of these two formulations will be provided in the paper accompanying fddm, but we will continue to use the traditional formulation for the remainder of this vignette.
Now allowing the drift rate to vary across trials (i.e. \(sv > 0\)), we should have two density functions. However, as only the small-time variable drift rate density function has been available in the literature (Blurton, Kesselmeier, and Gondan 2017), we provide the derivation of the large-time variable drift rate density function in the fddm paper. The large-time and small-time variable drift rate density function are: \[\begin{align}
f_\ell(t | v, \eta^2, a, w) &= \frac{\pi}{a^2 \sqrt{1 + \eta^2 t}}
\exp{ \left( \frac{\eta^2 a^2 w^2 -2vaw -v^2 t}{2 (1 + \eta^2 t)} \right)}
\sum_{j = 1}^{\infty} j \sin \left( j w \pi \right) \exp{ \left( -\frac{j^2 \pi^2 t}{2a^2} \right)}, \label{eq:var-l}\\
f_s(t | v, \eta^2, a, w) &= \frac{a}{\sqrt{2 \pi t^3 \left( 1 + \eta^2 t \right)}}
\exp{ \left( \frac{\eta^2 a^2 w^2 -2vaw -v^2 t}{2 (1 + \eta^2 t)} \right)}
\sum_{j = -\infty}^{\infty} (w + 2j) \exp{ \left( -\frac{a^2}{2t} \left( w + 2j \right)^2 \right)}. \label{eq:var-s}
\end{align}\]
Immediately of note is that the infinite summation for each time scale is the same regardless of the inclusion of variability in the drift rate. It then follows that there exists a term \(M\) such that the density function for the constant drift rate multiplied by \(M\) yields the density function for the variable drift rate. That is, \(M \cdot f(t | v, a, w) = f(t | v, a, w, \eta^2)\) from the above equations; this \(M\) works for converting both the large-time and small-time constant drift rate densities to variable drift rate densities. Although we do not use this term, it may be useful in adapting current algorithms to easily outputting the density with variable drift rate. Note that there are some issues with simply scaling the constant drift rate density, so please see the Validity Vignette for more information about the potential problems with this conversion. The multiplicative term \(M\) is given below: \[\begin{equation} M = \frac{1}{\sqrt{1 + \eta^2 t}} \exp \left( vaw + \frac{v^2 t}{2} + \frac{\eta^2 a^2 w^2 -2vaw -v^2 t}{2 (1 + \eta^2 t)} \right). \nonumber \end{equation}\]
The main issue of these families of density functions is that they all contain an infinite sum that must be approximated. Since there is no closed form analytical solution to this infinite sum, we instead calculate only a partial sum by truncating the sequence of terms after a certain point. We cannot actually calculate the true value of the density function, but we can mathematically prove that we can get arbitrarily close to the true value; the proof of this fact is provided in the paper accompanying the fddm package. The nature of this truncation has been the topic of many papers in the literature, but the underlying idea supporting all of the methods is the same: the user specifies an allowable error tolerance, and the algorithm calculates terms of the infinite sum until the result is within the allowed error tolerance of the true value.
The methods in the literature pre-calculate the number of terms required for the infinite sum to converge within the allowed error tolerance, and this number of terms is referred to as \(k_\ell\) and \(k_s\) for the large-time and small-time infinite sums, respectively. Navarro and Fuss (2009) include a method for calculating \(k_b\), the number of required terms for the infinite sum when combining the density functions of the two time scales. In addition to these existing methods, we add a new combination of large-time and small-time density functions and also provide a novel method that does not perform this pre-calculation. Note that in each method that pre-calculates the number of terms, the response time \(t\) is scaled inversely by \(a^2\), that is \(t' := \tfrac{t}{a^2}\). Also note that for the rest of this vignette, the ceiling function will be denoted by \(\lceil \cdot \rceil\).
The large-time density functions, Equations \(\eqref{eq:con-l}\) and \(\eqref{eq:var-l}\), have an infinite sum that runs for all of the positive integers. For a given error tolerance \(\epsilon\), Navarro and Fuss (2009) provide an expression for \(k_\ell\), the number of terms required for the large-time infinite sum to be within \(\epsilon\) of the true value of the density function. Thus the infinite sum becomes finite: \[\begin{equation} \label{eq:kl} \sum_{j = 1}^{k_\ell^\text{Nav}} j \sin \left( j w \pi \right) \exp{ \left( -\frac{j^2 \pi^2 t'}{2} \right)}. \nonumber \end{equation}\]
It remains to find the value of \(k_\ell^\text{Nav}\) that ensures the truncated sum is \(\epsilon\)-close to the true value. Navarro and Fuss (2009) provide a derivation in their paper that finds an upper bound for the tail of the sum, the sum of all terms greater than \(k_\ell^\text{Nav}\) (i.e., the error). Then they back-calculate the number of terms required to force this upper bound on the error to be less than \(\epsilon\), since then the actual error must also be less than \(\epsilon\). The resulting number of terms is: \[\begin{equation} \label{eq:kl-Nav} k_\ell^{\text{Nav}} \left( t', \epsilon \right) = \left\lceil \max \left\{ \sqrt{\frac{-2 \log(\pi t' \epsilon)}{\pi^2 t'}}, \frac{1}{\pi \sqrt{t'}} \right\} \right\rceil. \tag{L.1} \end{equation}\]
This method is often viewed as the most inefficient of the available options in the literature; however, this method proves to be extremely efficient in particular areas of the parameter space (typically for a large \(t'\)). The argument used to implement this method in dfddm is to set the parameter scale = "large" in the function call. In this case the other parameters, n_terms_small and summation_small, are ignored to avoid any potential clashes with overridden options.
The small-time approximations, Equations \(\eqref{eq:con-s}\) and \(\eqref{eq:var-s}\), also have an infinite sum, but this one runs over all of the integers – from negative infinity to positive infinity. Given this infinite nature in both directions, it is impossible to rigorously define the number of terms required to achieve the \(\epsilon\)-accuracy because we don’t know where to start counting the terms. To solve this issue, we rearrange the terms in the sum into the sequence \(\left\{ b_0, b_{-1}, b_1, \dots, b_{-j}, b_j, b_{-(j+1)}, b_{j+1}, \dots \right\}\); this allows us not only to count the terms in a sensible manner but also to define \(k_s\) as the index of the sequence where the truncation should occur. Then we can write the truncated version of the sum: \[\begin{equation} \label{eq:ks} \sum_{j = -k_s}^{k_s} (w + 2j) \exp{ \left( -\frac{a^2}{2t} \left( w + 2j \right)^2 \right)}. \nonumber \end{equation}\]
To choose the small-time methods when using dfddm, set the optional parameter scale = "small" in the function call. You can also set the optional parameter summation_small = "2017" or summation_small = "2014", but it is recommended to ignore this parameter so it retains its default value of “2017” that evaluates slightly faster than its counterpart. This parameter controls the style of summation used in the small-time approximation, and more details on the differences between these two styles can be found in the paper accompanying fddm. The final parameter, n_terms_small, will be discussed in the following three subsections.
Similarly to their large-time solution, Navarro and Fuss (2009) provide an expression for \(k_s\) given an error tolerance \(\epsilon\). They follow a similar idea to their large-time derivation by bounding the error of the truncation with integrals then back-calculating the number of terms required to keep the error less than \(\epsilon\). The resulting approximation yields: \[\begin{equation} \label{eq:ks-Nav} k_s^{\text{Nav}} \left( t', \epsilon \right) = \left\lceil \max \left\{ 2 + \sqrt{-2t' \log(2 \epsilon \sqrt{2 \pi t'})}, 1 + \sqrt{t'} \right\} \right\rceil. \tag{S.1} \end{equation}\]
To use this method, set n_terms_small = "Navarro" and keep scale = "small" in the function call. The parameter summation_small should be ignored so that it retains its default value to obtain the best performance.
Gondan, Blurton, and Kesselmeier (2014) improve on the approximation provided by Navarro and Fuss (2009) by using fewer terms in the sum to achieve the same level of accuracy, translating into a faster computation time. Given the error tolerance \(\epsilon\), Gondan, Blurton, and Kesselmeier (2014) exploit the intrinsic nature of the terms in the infinite sum to derive a smaller value for \(k_s\). After a few terms of the sum have been evaluated, each remaining term is smaller in absolute value than the previous term. They combine this knowledge with the observation that the terms alternate in sign to produce the following required number of terms: \[\begin{equation} \label{eq:ks-Gon} \begin{aligned} k_s^{\text{Gon}} \left( t', w, \epsilon \right) &= \left\lceil \max \left\{ \tfrac{1}{2} \left( \sqrt{2t'} - w \right), \tfrac{1}{2} \left( \sqrt{-t' (u_\epsilon - \sqrt{-2 u_\epsilon -2})} - w \right) \right\} \right\rceil,\\ u_\epsilon &= \min \left\{ -1, \log(2 \pi t'^2 \epsilon^2) \right\}. \end{aligned} \tag{S.2} \end{equation}\]
To use this method, set n_terms_small = "Gondan" and keep scale = "small" in the function call. The parameter summation_small should be ignored so that it retains its default value to obtain the best performance.
If we consider the terms of the infinite sum as the sequence defined above, the series alternates in sign \((+, -, +, ...)\); moreover, the series eventually decreases monotonically (in absolute value) due to the exponential term. Combining and exploiting these two mathematical properties has been the cornerstone of the previous approximations, but we will instead truncate the sum using a novel method. This method does not pre-calculate the number of terms required to achieve the given error tolerance. Instead, the general idea of this method is to take full advantage of the alternating and decreasing nature of the terms in the infinite sum by applying a handy theorem (commonly known as Leibniz’s rule) to place an upper bound on the truncation error after including so many terms. It has been proven that this upper bound is in fact the absolute value of the next term in the sequence, thus we can truncate the infinite sum once one of its terms is less than the desired error tolerance, \(\epsilon\). Hence we do not consider the number of terms in the sum, rather just that the terms in the summation will eventually be small enough. The validity of this method is proven in the paper that accompanies the fddm package.
To use this method, set n_terms_small = "SWSE" and keep scale = "small" in the function call. The parameter summation_small should be ignored so that it retains its default value to obtain the best performance.
A sensible next approach to approximating the density of the DDM is to use some combination of the large-time and small-time density functions. As their names suggest, each density function approximation has a section of the parameter space where it outperforms the other one. Essentially these methods involve calculating the number of terms required for both the large-time and small-time density functions, then using whichever approximation requires fewer terms. The goal is to use each approximation where it is efficient and avoid the areas of the parameter space where it is horribly inefficient.
To use this method of evaluating the PDF, set the optional parameter scale = "both" in the function call. These methods include similar options to the small-time approximations from the previous section, and thus it is possible to set the optional parameters in a similar fashion to above. The value of summation_small triggers the same effect as above since it only affects the small-time approximation to the PDF; we recommend ignoring this parameter so it assumes its default value of “2017” that evaluates slightly faster than its counterpart. The final parameter, n_terms_small, will be discussed in the following two subsections. Note that there is no “SWSE” option for scale = "both", and setting n_terms_small = "SWSE" will override the scale = "both" setting. As there is only one option for the large time approximation, there are no optional parameters to set for this part of the combined time scale approximation.
Navarro and Fuss (2009) initially suggested a method where you calculate both \(k_s^\text{Nav}\) and \(k_l^\text{Nav}\) and use whichever one has a smaller value. However, one issue with their original method arises when the user inputs a vector of response times. They calculate the maximum value of \(k_s^\text{Nav}\) for all of the input response times and then compare that to the maximum value of \(k_l^\text{Nav}\) for all of the input response times. They end up with one \(k_b^\text{Nav}\) for all of the input response times, and this one value can be higher than normal for a response time due to the maximization of the approximation over all of the response times. This over-accuracy is not necessarily a bad thing, but it does mean that it is possible to get different densities for the same input response time depending on the other response times that have been input. To correct for this, we simply do not maximize \(k_s^\text{Nav}\) or \(k_l^\text{Nav}\) over all of the input response times and instead calculate one value of \(k_b^\text{Nav}\) for each input response time. \[\begin{equation} \label{eq:kb-Nav} k_b^\text{Nav} \left( t', w, \epsilon \right) = \min \left\{ k_s^{\text{Nav}}, k_l^{\text{Nav}} \right\}. \tag{B.1} \end{equation}\]
To use this method, set n_terms_small = "Navarro" and keep scale = "both" in the function call. The parameter summation_small should be ignored so that it retains its default value to obtain the best performance.
This combination of methods has not been explored in the literature before, but it works very similarly to the combination above. The only difference is that we use the Gondan, Blurton, and Kesselmeier (2014) approximation for the small-time instead of the one provided by Navarro and Fuss (2009). Since \(k_s^\text{Gon} \leq k_s^\text{Nav}\), this method will give approximations at least as fast as the above combination. \[\begin{equation} \label{eq:kb-Gon} k_b^{\text{Gon}} \left( t', w, \epsilon \right) = \min \left\{ k_s^{\text{Gon}}, k_l^{\text{Nav}} \right\}. \tag{B.2} \end{equation}\]
To use this method, set n_terms_small = "Gondan" and keep scale = "both" in the function call. The parameter summation_small should be ignored so that it retains its default value to obtain the best performance.
sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18363)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=English_United Kingdom.1252
#> [3] LC_MONETARY=English_United Kingdom.1252
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United Kingdom.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplot2_3.3.2.9000 reshape2_1.4.4 microbenchmark_1.4-7
#> [4] RWiener_1.3-3 rtdists_0.11-2 fddm_0.1-1
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.4.6 pillar_1.4.4 compiler_4.0.2 plyr_1.8.6
#> [5] tools_4.0.2 digest_0.6.25 evd_2.3-3 evaluate_0.14
#> [9] lifecycle_0.2.0 tibble_3.0.1 gtable_0.3.0 lattice_0.20-41
#> [13] pkgconfig_2.0.3 rlang_0.4.6 Matrix_1.2-18 yaml_2.2.1
#> [17] mvtnorm_1.1-1 expm_0.999-4 xfun_0.15 withr_2.2.0
#> [21] dplyr_1.0.0 stringr_1.4.0 knitr_1.29 generics_0.0.2
#> [25] vctrs_0.3.1 tidyselect_1.1.0 grid_4.0.2 ggnewscale_0.4.1
#> [29] glue_1.4.1 R6_2.4.1 survival_3.2-3 rmarkdown_2.3
#> [33] farver_2.0.3 purrr_0.3.4 magrittr_1.5 scales_1.1.1
#> [37] htmltools_0.5.0 ellipsis_0.3.1 splines_4.0.2 colorspace_1.4-1
#> [41] labeling_0.3 stringi_1.4.6 gsl_2.1-6 munsell_0.5.0
#> [45] msm_1.6.8 crayon_1.3.4Blurton, Steven P, Miriam Kesselmeier, and Matthias Gondan. 2017. “The First-Passage Time Distribution for the Diffusion Model with Variable Drift.” Journal of Mathematical Psychology 76: 7–12.
Feller, William. 1968. An Introduction to Probability Theory and Its Applications: Volume I. Vol. 1. John Wiley & Sons.
Gondan, Matthias, Steven P Blurton, and Miriam Kesselmeier. 2014. “Even Faster and Even More Accurate First-Passage Time Densities and Distributions for the Wiener Diffusion Model.” Journal of Mathematical Psychology 60: 20–22.