dfddmFunction 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. 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),
\(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 are several different methods for approximating the PDF of the
DDM, and there are three optional control parameters in
dfddm() that can be used to indicate which method should be
used in the function call: switch_mech,
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 implementation 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 and most stable algorithm. Note that
precedence for the optional parameters is first given to checking if the
default implementation is
selected. If not, precedence is then given to the
switch_mech parameter value; for example,
switch_mech = "large" will ignore the
summation_small input value.
Note that there are actually two probability density functions of the DDM: the PDF with respect to the upper boundary, and the PDF with respect to the lower boundary. Following the precedent set by the literature, we use the general terminology “PDF” or “density function” to mean the probability density function with respect to the lower boundary. Should the probability density 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 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 by truncating the infinite sum. This vignette details the various methods used in the literature to approximate the infinite sum in the density function for the DDM.
The author of 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. Each of these forms contains an infinite sum, and they are known individually 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. The improved speed is due to the ease at which the infinite sums can be adequately approximated with a finite truncation. In other words, the efficiency of each approximation is dictated by the number of terms required in the truncated sum to achieve the pre-specified precision; fewer required terms translates to fewer computations and thus generally faster computation time.
When the drift rate is held constant (i.e. \(\eta = 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 are 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. \(\eta > 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 et al. 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, 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, 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 value \(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 so that the truncated 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_c\), 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 novel method that does not perform this pre-calculation, and we also provide two new combinations of the large-time and small-time density functions. 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
large \(t'\)). To implement this
method in dfddm() (and completely excluding the small-time
methods), the user must set the parameter
switch_mech = "large" in the function call; in this case,
the other parameters are ignored.
The small-time approximations, Equations \(\eqref{eq:con-s}\) and \(\eqref{eq:var-s}\), also contain an infinite sum, but this sum 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 = -\left\lfloor\frac{k_s}{2}\right\rfloor}^{\left\lfloor\frac{k_s}{2}\right\rfloor} (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() (and
completely excluding the large-time method), set the optional parameter
switch_mech = "small" in the function call. You can also
set the optional control 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 control parameter,
n_terms_small, will be discussed in the following three
subsections.
After Navarro and Fuss (2009) published their paper, Gondan et al. (2014) introduced another method for calculating the required number of terms in the truncated small-time summation. It is important to note, however, that Gondan et al. (2014) provided the number of required pairs of terms in the \(S_{14}\) summation style, and not the number of required individual terms. As we want the number of individual terms, we adapt their formula and define it given a desired precision \(\epsilon\): \[\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 switch_mech = "small" and
n_terms_small = "Gondan" 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 \((+, -,
+, \dots)\); 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 method suggested by Gondan et al.
(2014). 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 the alternating series test) 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 (in absolute value) 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 swich_mech = "small" and
n_terms_small = "SWSE" 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 the approximations perform poorly. The main
control parameter used to indicate this preference is
switch_mech.
Fundamentally, each combined time scale method uses either the
large-time or small-time approximation for each calculation of the
density function. As there is only one option for the large-time
approximation, there are no optional control parameters to set for this
part of the combined time scale approximation. In contrast, there are
multiple options for the small-time part of the combined time scale
approximation. However, the effect of the control parameter
summation_small is consistent throughout the small-time
methods; we recommend leaving it to its default value of “2017” for the
best performance. The remaining subsections of this vignette detail how
to set the rest of the control parameters in the function call to use a
particular set of methods for calculating the combined time scale
density function.
This combination of methods has not been explored in the literature before, but it works very similarly to the Navarro-Navarro combination above. The only difference is that we use the Gondan et al. (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 should never be less efficient than the previous combined time scale method. \[\begin{equation} \label{eq:kc-Gon} k_c^{\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 switch_mech = "terms" and
n_terms_small = "Gondan" in the function call. The
parameter summation_small should be ignored so that it
retains its default value to obtain the best performance.
The SWSE approximation to the small-time density function differs
from the Navarro or Gondan approximations in that it does not
pre-calculate \(k_s\), the number of
terms in the infinite sum that are required to achieve the desired
precision; instead, the infinite sum is truncated when the individual
terms of the sum become “small enough.” This method of truncating the
infinite sum poses a problem of how to incorporate this method with the
Navarro large-time approximation that relies on pre-calculating \(k_\ell^{\text{Nav}}\). We will introduce
two heuristics for determining when to use the small-time approximation
or the large-time approximation. Both heuristics will use a fourth
parameter called switch_thresh, but each heuristic will use
the new parameter in a slightly different way. The validity of these two
methods is proven in the paper that accompanies the fddm
package. Note that in this paper, the parameter
max_terms_large is labelled \(\delta\).
The first heuristic compares switch_thresh to \(k_\ell^{\text{Nav}}\); in this case,
switch_thresh is treated as the proxy for the required
number of terms for approximating the truncated small-time infinite sum.
If \(k_\ell^{\text{Nav}}\) \(\le\) switch_thresh, then the
Navarro large-time approximation is used. On the other hand, if \(k_\ell^{\text{Nav}}\) \(>\) switch_thresh, then the
SWSE small-time approximation is used. This method essentially checks
the efficiency of the large-time approximation relative to
switch_thresh. The user can alter the behavior of this
method by setting the optional parameter switch_thresh to
any non-negative integer; the default value for this parameter when
using this method is \(1\).
To use this method, set switch_mech = "terms_large". The
user may wish to specify a particular threshold (measured in the
required number of terms for the truncated sum) for switching between
the small-time and large-time approximations by setting the parameter
switch_thresh (i.e., switch_thresh = 1). The
parameter summation_small should be ignored so that it
retains its default value to obtain the best performance.
The second heuristic avoids all comparison with \(k_\ell^{\text{Nav}}\) as it simply
considers the effective response time, \(t' := \tfrac{t}{a^2}\), to be the
indicator of a small or large response time. \(t'\) is compared to a new parameter
switch_thresh so that if \(t'\) \(>\) switch_thresh, then the
Navarro large-time approximation is used; otherwise, the SWSE small-time
approximation is used. The user can alter the behavior of this method by
setting the optional parameter switch_thresh to any
non-negative real number; the default value for this parameter when
using this method is \(0.8\).
Since this is the default implementation, all four of the optional
parameters (switch_mech, switch_thresh,
n_terms_small, and summation_small) can be
ignored. The user may wish to specify a particular threshold (measured
in seconds) for switching between the small-time and large-time
approximations by setting the parameter switch_thresh
(i.e., switch_thresh = 0.8). We recommend ignoring the
parameter summation_small so that it retains its default
value and achieves optimal performance. Note that the parameter
n_terms_small is ignored because the only option is “SWSE”.
sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggforce_0.5.0 ggplot2_4.0.3 reshape2_1.4.5
#> [4] microbenchmark_1.5.0 RWiener_1.3-3 rtdists_0.11-5
#> [7] fddm_1.0-2 rmarkdown_2.31
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 generics_0.1.4 stringi_1.8.7 lattice_0.22-9
#> [5] digest_0.6.39 magrittr_2.0.5 estimability_1.5.1 evaluate_1.0.5
#> [9] grid_4.6.0 RColorBrewer_1.1-3 mvtnorm_1.3-7 fastmap_1.2.0
#> [13] plyr_1.8.9 jsonlite_2.0.0 Matrix_1.7-5 ggnewscale_0.5.2
#> [17] Formula_1.2-5 survival_3.8-6 scales_1.4.0 tweenr_2.0.3
#> [21] codetools_0.2-20 jquerylib_0.1.4 cli_3.6.6 rlang_1.2.0
#> [25] expm_1.0-0 polyclip_1.10-7 gsl_2.1-9 splines_4.6.0
#> [29] withr_3.0.2 cachem_1.1.0 yaml_2.3.12 tools_4.6.0
#> [33] buildtools_1.0.0 vctrs_0.7.3 R6_2.6.1 emmeans_2.0.3
#> [37] lifecycle_1.0.5 stringr_1.6.0 MASS_7.3-65 pkgconfig_2.0.3
#> [41] pillar_1.11.1 bslib_0.11.0 gtable_0.3.6 glue_1.8.1
#> [45] Rcpp_1.1.1-1.1 tidyselect_1.2.1 xfun_0.57 tibble_3.3.1
#> [49] sys_3.4.3 knitr_1.51 msm_1.8.2 farver_2.1.2
#> [53] htmltools_0.5.9 labeling_0.4.3 maketools_1.3.2 compiler_4.6.0
#> [57] S7_0.2.2 evd_2.3-7.1