Description of Methods in pfddm

Kendal Foster and Henrik Singmann

July 01, 2024

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.

Mathematical Background

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.

The Distribution Function Approximations

There are two main ways of calculating the CDF of the DDM: approximating the analytic form of CDF itself (Blurton, Kesselmeier, and Gondan 2017), and numerically solving a partial differential equation (PDE) whose solution is the CDF (Voss and Voss 2008). 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.

Infinite Sum Methods

Similarly to the density function of the DDM, there is a large-time variant and a small-time variant of the CDF for a constant drift rate (i.e., \(\eta = 0\)) (see Blurton, Kesselmeier, and Gondan (2012) for details). There is no analytic large-time variant for \(\eta > 0\) because the resulting integral is incalculable (see Tuerlinckx (2004) for more details). However, Blurton, Kesselmeier, and Gondan (2017) provide the small-time variant of the CDF for \(\eta > 0\); this variant contains an infinite sum that must be approximated. Blurton, Kesselmeier, and Gondan (2017) even provide two different forms of the small-time variant of the CDF: one that uses the standard normal CDF in each term of the infinite sum, and one that instead utilizes Mills ratio in each term of the infinite sum.

Similarly to the small-time variant of the PDF, the small-time variant of the CDF contains an infinite sum whose terms alternate and decay to zero. Like in the SWSE method of dfddm(), we can exploit these mathematical properties and bound the approximation error. In this case, the approximation error, \(\epsilon\), is bounded above by the absolute value of the next term in the sum. For example, if the sum includes \(k\) terms (\(b_1 + b_2 + \dots + b_k\)), then the approximation error is at most the absolute value of the \(k+1\)th term (\(\epsilon \leq \left\lvert b_{k+1} \right\rvert\)). These mathematical properties allow us to explicitly pre-define and control the precision of the calculation before actually doing it.

Moreover, in the Benchmark Results below, it is shown that this class of approximation methods is quite fast. The combination of directly controlled precision and fast calculation times are why pfddm() uses the infinite sum methods instead of the PDE method for approximating the CDF of the DDM.

The following two subsections detail the two slightly different variations on the small-time CDF provided by Blurton, Kesselmeier, and Gondan (2017).

Normal CDF

When Blurton, Kesselmeier, and Gondan (2017) derived the the small-time CDF with inter-trial variability in the drift rate (i.e., \(\eta > 0\)), they first arrived at the following equation for the distribution function1:

\[\begin{align} \label{eq:cdf_ncdf} F(t ~|~ v, \eta, a, w) = &\exp{\left( \frac{1}{2} \eta^2 a^2 w^2 - v a w \right)} \times\\ &\sum_{j = 0}^{\infty} (-1)^j \exp{\left( \frac{1}{2} \eta^2 r_j^2 \right)} \Bigg( \exp{(-\gamma r_j)} \Phi \left(\frac{\gamma t - \lambda r_j}{\rho} \right) + \exp{(\gamma r_j)} \Phi \left(\frac{-\gamma t - \lambda r_j}{\rho} \right) \Bigg)\nonumber \end{align}\]

where

\[\begin{align*} r_j &= \begin{cases} a \big( j + w \big) & \text{if } j \text{ is even,}\\ a \big( j + (1-w) \big) & \text{if } j \text{ is odd,} \end{cases}\\[0.5ex] \gamma &= v - \eta^2 a w,\\[0.1ex] \lambda &= 1 + \eta^2 t,\\[0.1ex] \rho &= \sqrt{\lambda t}, \nonumber \end{align*}\]

and \(\Phi(x)\) is the standard normal CDF evaluated at \(x\).

One important note is that the \(\exp{\left( \frac{1}{2} \eta^2 r_j^2 \right)}\) expression in each term of the infinite sum tends to infinity as \(j \to \infty\). This is balanced by the terms inside the parentheses decaying to \(0\) quickly, yet it is still cause for concern from a computational perspective. If the growing exponential overflows to positive infinity before the sum is truncated, then the approximation is no longer useful2.

To use this method, set method = "NCDF" in the function call to pfddm(). We caution use of this approximation method, as the default-method (using Mills ratio) is both faster and more stable.

Mills Ratio

After deriving Equation \(\ref{eq:cdf_ncdf}\), Blurton, Kesselmeier, and Gondan (2017) manipulated their result to use the Mills ratio for its favorable numerical properties (see Blurton, Kesselmeier, and Gondan (2012) for more discussion on these properties of Mills ratio). The Mills ratio is defined as

\[\begin{equation} \label{eq:mills} M[x] = \frac{1 - \Phi(x)}{\phi(x)} = \frac{\Phi(-x)}{\phi(x)}, \end{equation}\]

where \(\Phi(x)\) is the standard normal CDF evaluated at \(x\), and \(\phi(x)\) is the standard normal PDF evaluated at \(x\). Then the CDF of the DDM becomes

\[\begin{align} \label{eq:cdf_mills} F(t ~|~ v, \eta, a, w) = &\exp{ \left( \frac{\eta^2 a^2 w^2 - 2 v a w - v^2 t}{2 \lambda} \right) } \times\\ &\sum_{j = 0}^{\infty} (-1)^j \phi \left( \frac{r_j}{\sqrt{t}} \right) \Bigg( M \left( \frac{\lambda r_j - \gamma t}{\rho} \right) + M \left( \frac{\lambda r_j + \gamma t}{\rho} \right) \Bigg)\nonumber \end{align}\]

where again

\[\begin{align*} r_j &= \begin{cases} a \big( j + w \big) & \text{if } j \text{ is even,}\\ a \big( j + (1-w) \big) & \text{if } j \text{ is odd,} \end{cases}\\[0.5ex] \gamma &= v - \eta^2 a w,\\[0.1ex] \lambda &= 1 + \eta^2 t,\\[0.1ex] \rho &= \sqrt{\lambda t}. \nonumber \end{align*}\]

Whereas this form does not have the issue of the exponential increasing to infinity, the Mills ratio does introduce a potential divide-by-zero error if one is not careful with handling large arguments in the PDF in the denominator of the Mills ratio3.

Fortunately, there exists a useful approximation to the Mills ratio for large arguments. Abramowitz and Stegun (1964) provide an asymptotic expansion for the complementary CDF of the standard normal distribution in Expression 26.2.13 in their book; dividing this asymptotic expansion by the PDF of the standard normal distribution yields the Mills ratio. This approximation was found in the zeta function in the sn package (Azzalini 2021) for R, which in turn was recommended by Gondan, Blurton, and Kesselmeier (2014). The zeta function actually calculates the inverse Mills ratio, so we use the reciprocal of the approximation used in the zeta function.

By combining the standard Mills ratio with this approximation, we can achieve a robust algorithm that can be used to consistently calculate the CDF. Moreover, this approximation method is quite fast, as demonstrated in the Benchmark Results below. The speed and stability of this algorithm earns our recommendation and is the default method in pfddm().

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

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 Voss and Voss (2008).

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, 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 Voss and Voss (2008) (using the function rtdists::pdiffusion() from the 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 Voss and Voss (2008) 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 (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 package, as this function uses the PDE method to calculate the CDF.

Benchmarking the Distribution Function Approximations

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 Blurton, Kesselmeier, and Gondan (2017), 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")

Vectorized Benchmark Data and Results

The first benchmark test will record the algorithms’ performances across the parameter space, with the response times input as a vector. Inputting the response times in this way is both more common (as it is typically easier to input a vector into a function as opposed to writing a loop) and beneficial for the R-based implementation from Blurton, Kesselmeier, and Gondan (2017) to exploit R’s vectorization (making this implementation as fast as possible).

We write the function that we will use to systematically benchmark each implementation at each combination of parameters.

rt_benchmark_vec <- function(RT, resp, V, A, t0 = 1e-4, W = 0.5, SV = 0.0,
                             err_tol = 1e-6, times = 1000) {

  fnames <- c("Mills", "NCDF", "Blurton", "rtdists")
  nf <- length(fnames) # number of functions being benchmarked
  nV <- length(V)
  nA <- length(A)
  nW <- length(W)
  nSV <- length(SV)
  resp <- rep(resp, length(RT)) # for RWiener

  # Initialize the dataframe to contain the microbenchmark results
  mbm_res <- data.frame(matrix(ncol = 4+nf, nrow = nV*nA*nW*nSV))
  colnames(mbm_res) <- c('V', 'A', 'W', 'SV', fnames)
  row_idx <- 1

  # Loop through each combination of parameters and record microbenchmark results
  for (v in 1:nV) {
    for (a in 1:nA) {
      for (w in 1:nW) {
        for (sv in 1:nSV) {
          mbm <- microbenchmark(
          Mills = pfddm(rt = RT, response = resp, a = A[a], v = V[v], t0 = t0,
                        w = W[w], sv = SV[sv], log = FALSE, method = "Mills",
                        err_tol = err_tol),
          NCDF = pfddm(rt = RT, response = resp, a = A[a], v = V[v], t0 = t0,
                       w = W[w], sv = SV[sv], log = FALSE, method = "NCDF",
                       err_tol = err_tol),
          Blurton = G_0(t = RT-t0, a = A[a], nu = V[v], w = W[w],
                        eta2 = SV[sv]*SV[sv], sigma2 = 1, eps = err_tol), # only "lower" resp
          rtdists = pdiffusion(RT, resp, a = A[a], v = V[v], t0 = t0,
                               z = W[w]*A[a], sv = SV[sv], precision = 3),
          times = times)
        # add the v, a, and w values to the dataframe
        mbm_res[row_idx, 1] <- V[v]
        mbm_res[row_idx, 2] <- A[a]
        mbm_res[row_idx, 3] <- W[w]
        mbm_res[row_idx, 4] <- SV[sv]
        # add the median microbenchmark results to the dataframe
        for (i in 1:nf) {
          mbm_res[row_idx, 4+i] <- median(mbm[mbm[,1] == fnames[i],2])
        }
        # iterate start value
        row_idx <- row_idx + 1
        }
      }
    }
  }
  return(mbm_res)
}

We will not actually run the benchmark tests in this vignette as it can take a while to run, and instead we will load pre-run benchmark data before we plot the results.

bm_vec <- rt_benchmark_vec(RT = RT, resp = "lower", V = V, A = A, t0 = t0,
                           W = W, SV = SV, err_tol = err_tol, times = 1000)

Here, we load the pre-run benchmark data and plot the results as side-by-side violin plots. The horizontal axis shows the implementation, and the vertical axis shows the benchmark time. Each individual violin plot shows a mirrored density estimate; overlaying the violin plot, the boxplot shows the median in addition to the first and third quartiles; the horizontal dashed line shows the mean.

# load data, will be in the variable 'bm_vec'
load(system.file("extdata", "pfddm_distribution", "bm_vec_0-2.Rds",
                 package = "fddm", mustWork = TRUE))

t_idx <- match("SV", colnames(bm_vec))
bm_vec[, -seq_len(t_idx)] <- bm_vec[, -seq_len(t_idx)]/1000 # convert to microseconds
mbm_vec <- melt(bm_vec, measure.vars = -seq_len(t_idx),
                variable.name = "FuncName", value.name = "time")

Names_vec <- c("Mills", "NCDF", "Blurton", "rtdists")
Color_vec <- c("#b34d4d", "#4d80b3", "#c5a687", "#ac8053")
Outline_vec <- c("#b34d4d", "#4d80b3", "#c5a687", "#ac8053")

mi <- min(bm_vec[, (t_idx+1):(ncol(bm_vec)-2)])
ma <- max(bm_vec[, (t_idx+1):(ncol(bm_vec)-2)])

ggplot(mbm_vec, aes(x = factor(FuncName, levels = Names_vec), y = time,
                    color = factor(FuncName, levels = Names_vec),
                    fill = factor(FuncName, levels = Names_vec))) +
  geom_violin(trim = TRUE, alpha = 0.5) +
  scale_color_manual(values = Outline_vec, guide = FALSE) +
  scale_fill_manual(values = Color_vec, guide = FALSE) +
  geom_boxplot(width = 0.15, fill = "white", alpha = 0.5) +
  stat_summary(fun = mean, geom = "errorbar",
               aes(ymax = ..y.., ymin = ..y..),
               width = .35, linetype = "dashed") +
  scale_x_discrete(labels = c(
    bquote(F[s] ~ Mills), bquote(F[s] ~ NCDF), "Blurton (Mills)", "rtdists")) +
  facet_zoom(ylim = c(mi, ma)) +
  labs(x = "Implementation", y = "Time (microseconds)") +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16, angle = 90,
                                   vjust = 0.5, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20,
                                    margin = margin(10, 0, 0, 0)),
        axis.title.y = element_text(size = 20,
                                    margin = margin(0, 10, 0, 0)),
        legend.position = "none")
#> Warning: The dot-dot notation (`..y..`) was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(y)` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
#> ggplot2 3.3.4.
#> ℹ Please use "none" instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

Depending on your system, the results might differ from the ones shown here, but they should replicate the general pattern. The two implementations in pfddm() are both on average faster than the pure R code implementation and the implementation in rtdists, but each implementation still has quite a long tail on its distribution of benchmark times. These inefficiencies will be explored in the next section, where we can isolate the response time inputs and determine in which parts of the parameter space these algorithms struggle. Of the two implementations in pfddm(), the one utilizing the Mills ratio is on average faster than the one that just uses the standard normal CDF. Although the Mills ratio must calculate an additional standard normal PDF (i.e., the denominator of the Mills ratio), the useful approximation from Abramowitz and Stegun (1964) allows for a fast and accurate approximation to it.

Non-vectorized Benchmark Data and Results

The second benchmark test will record the algorithms’ performances across the parameter space, with the response times input individually (i.e., not as a vector). Inputting the response times in this way are both more common (as it is typically easier to input a vector into a function as opposed to writing a loop) and beneficial for the R-based implementation from Blurton, Kesselmeier, and Gondan (2017) to exploit R’s vectorization (making this implementation as fast as possible).

We write the function that we will use to systematically benchmark each implementation at each combination of parameters.

rt_benchmark_ind <- function(RT, resp, V, A, t0 = 1e-4, W = 0.5, SV = 0.0,
                             err_tol = 1e-6, times = 100) {
  fnames <- c("Mills", "NCDF", "Blurton", "rtdists")
  nf <- length(fnames) # number of functions being benchmarked
  nRT <- length(RT)
  nV <- length(V)
  nA <- length(A)
  nW <- length(W)
  nSV <- length(SV)

  # Initialize the dataframe to contain the microbenchmark results
  mbm_res <- data.frame(matrix(ncol = 5+nf, nrow = nRT*nV*nA*nW*nSV))
  colnames(mbm_res) <- c('RT', 'V', 'A', 'W', 'SV', fnames)
  row_idx <- 1

  # Loop through each combination of parameters and record microbenchmark results
  for (rt in 1:nRT) {
    for (v in 1:nV) {
      for (a in 1:nA) {
        for (w in 1:nW) {
          for (sv in 1:nSV) {
            mbm <- microbenchmark(
            Mills = pfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], sv = SV[sv], log = FALSE,
                          method = "Mills", err_tol = err_tol),
            NCDF = pfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                         t0 = t0, w = W[w], sv = SV[sv], log = FALSE,
                         method = "NCDF", err_tol = err_tol),
            Blurton = G_0(t = RT[rt]-t0, a = A[a], nu = V[v], w = W[w],
                          eta2 = SV[sv]*SV[sv], sigma2 = 1, eps = err_tol), # only "lower" resp
            rtdists = pdiffusion(RT[rt], resp, a = A[a], v = V[v], t0 = t0,
                                 z = W[w]*A[a], sv = SV[sv], precision = 3),
            times = times)
          # add the v, a, and w values to the dataframe
          mbm_res[row_idx, 1] <- RT[rt]
          mbm_res[row_idx, 2] <- V[v]
          mbm_res[row_idx, 3] <- A[a]
          mbm_res[row_idx, 4] <- W[w]
          mbm_res[row_idx, 5] <- SV[sv]
          # add the median microbenchmark results to the dataframe
          for (i in 1:nf) {
            mbm_res[row_idx, 5+i] <- median(mbm[mbm[,1] == fnames[i],2])
          }
          # iterate start value
          row_idx <- row_idx + 1
          }
        }
      }
    }
  }
  return(mbm_res)
}

We will not actually run the benchmark tests in this vignette as it can take a while to run, and instead we will load pre-run benchmark data before we plot the results.

bm_ind <- rt_benchmark_ind(RT = RT, resp = "lower", V = V, A = A, t0 = t0,
                           W = W, SV = SV, err_tol = err_tol, times = 1000)

Here, we load the pre-run benchmark data and plot the results as a series of plots that show how the benchmark times change as a result of varying the response time input to the CDF approximation. On the horizontal axis is the response time (input to the CDF approximation), and the vertical axis displays the benchmark time; note that the vertical axis is not fixed across panels. The dark line in each panel indicates the mean benchmark time; the darker shaded region in each panel shows the 10% and 90% quantiles; and the lightly shaded region shows the minimum and maximum benchmark times.

# load data, will be in the variable 'bm_ind'
load(system.file("extdata", "pfddm_distribution", "bm_ind_0-2.Rds",
                 package = "fddm", mustWork = TRUE))

t_idx <- match("SV", colnames(bm_ind))
bm_ind[,-seq_len(t_idx)] <- bm_ind[, -seq_len(t_idx)]/1000 # convert to microseconds
mbm_ind <- melt(bm_ind, measure.vars = -seq_len(t_idx),
                variable.name = "FuncName", value.name = "time")

Names_meq <- c("Mills", "NCDF", "Blurton", "rtdists")
Color_meq <- c("#b34d4d", "#4d80b3", "#c5a687", "#ac8053")

my_labeller <- as_labeller(c(Mills = "F[s] ~ Mills",
                             NCDF = "F[s] ~ NCDF",
                             fl_Nav_09 = "f[l] ~ Nav",
                             Blurton = "Blurton (Mills)",
                             rtdists = "rtdists"),
                           default = label_parsed)

ggplot(mbm_ind, aes(x = RT, y = time,
                    color = factor(FuncName, levels = Names_meq),
                    fill = factor(FuncName, levels = Names_meq))) +
  stat_summary(fun.min = min, fun.max = max,
               geom = "ribbon", color = NA, alpha = 0.1) +
  stat_summary(fun.min = function(z) { quantile(z, 0.1) },
               fun.max = function(z) { quantile(z, 0.9) },
               geom = "ribbon", color = NA, alpha = 0.2) +
  stat_summary(fun = mean, geom = "line") +
  scale_x_log10(breaks = c(0.1, 0.25, 0.5, 1, 2),
                labels = as.character(c(0.1, 0.25, 0.5, 1, 2))) +
  scale_color_manual(values = Color_meq) +
  scale_fill_manual(values = Color_meq) +
  labs(subtitle = paste(
         "The darker shaded regions represent the 10% and 90% quantiles",
         "The lighter shaded regions represent the min and max times",
         sep = ";\n"),
       x = bquote(t ~ ", response time"),
       y = "Time (microseconds)") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.subtitle = element_text(size = 16,
                                     margin = margin(0, 0, 15, 0)),
        axis.text.x = element_text(size = 16, angle = 90,
                                   vjust = 0.5, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20,
                                    margin = margin(10, 0, 0, 0)),
        axis.title.y = element_text(size = 20,
                                    margin = margin(0, 10, 0, 0)),
        strip.text = element_text(size = 16),
        strip.background = element_rect(fill = "white"),
        legend.position = "none") +
  facet_wrap(~ factor(FuncName, levels = Names_meq), scales = "free_y",
             labeller = my_labeller)

This series of plots confirms that the inefficiencies of these algorithms occur when large response times are input, although the inefficiencies are not as dramatic as in the density function benchmark results. As there is no analytic form for the large-time distribution function that includes inter-trial variability in the drift rate, we must make do with only the small-time distribution function or the PDE method. As we discussed earlier in this vignette, approximating the infinite sum in the small-time distribution function is preferable to using the PDE method because it allows us to directly control the approximation error. Moreover, the small-time implementations in pfddm() are on average faster than the implementation of the PDE method in rtdists; they are also faster than the implementation in the pure R code from Blurton, Kesselmeier, and Gondan (2017).

R Session Info

sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Linux Mint 21.3
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
#>  [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: Europe/London
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggforce_0.4.2         ggplot2_3.5.1         reshape2_1.4.4       
#> [4] microbenchmark_1.4.10 RWiener_1.3-3         rtdists_0.11-5       
#> [7] fddm_1.0-2           
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.9         utf8_1.2.4         generics_0.1.3     stringi_1.8.4     
#>  [5] lattice_0.22-6     digest_0.6.35      magrittr_2.0.3     estimability_1.5.1
#>  [9] evaluate_0.24.0    grid_4.4.1         mvtnorm_1.2-5      fastmap_1.2.0     
#> [13] plyr_1.8.9         jsonlite_1.8.8     Matrix_1.7-0       ggnewscale_0.4.10 
#> [17] Formula_1.2-5      survival_3.7-0     fansi_1.0.6        scales_1.3.0      
#> [21] tweenr_2.0.3       codetools_0.2-20   jquerylib_0.1.4    cli_3.6.2         
#> [25] crayon_1.5.2       rlang_1.1.4        expm_0.999-9       polyclip_1.10-6   
#> [29] gsl_2.1-8          munsell_0.5.1      splines_4.4.1      withr_3.0.0       
#> [33] cachem_1.1.0       yaml_2.3.8         tools_4.4.1        dplyr_1.1.4       
#> [37] colorspace_2.1-0   vctrs_0.6.5        R6_2.5.1           emmeans_1.10.2    
#> [41] lifecycle_1.0.4    stringr_1.5.1      MASS_7.3-61        pkgconfig_2.0.3   
#> [45] pillar_1.9.0       bslib_0.7.0        gtable_0.3.5       glue_1.7.0        
#> [49] Rcpp_1.0.12        highr_0.11         xfun_0.45          tibble_3.2.1      
#> [53] tidyselect_1.2.1   knitr_1.47         msm_1.7.1          xtable_1.8-4      
#> [57] farver_2.1.2       htmltools_0.5.8.1  labeling_0.4.3     rmarkdown_2.27    
#> [61] compiler_4.4.1     evd_2.3-7

References

Abramowitz, Milton, and Irene A Stegun. 1964. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Vol. 55. US Government printing office.

Azzalini, A. 2021. The R Package sn: The Skew-Normal and Related Distributions Such as the Skew-\(t\) and the Sun (Version 2.0.0). Università di Padova, Italia. http://azzalini.stat.unipd.it/SN/.

Blurton, Steven P, Miriam Kesselmeier, and Matthias Gondan. 2012. “Fast and Accurate Calculations for Cumulative First-Passage Time Distributions in Wiener Diffusion Models.” Journal of Mathematical Psychology 56 (6): 470–75.

———. 2017. “The First-Passage Time Distribution for the Diffusion Model with Variable Drift.” Journal of Mathematical Psychology 76: 7–12.

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.

Tuerlinckx, Francis. 2004. “The Efficient Computation of the Cumulative Distribution and Probability Density Functions in the Diffusion Model.” Behavior Research Methods, Instruments, & Computers 36 (4): 702–16.

Voss, Andreas, and Jochen Voss. 2008. “A Fast Numerical Algorithm for the Estimation of Diffusion Model Parameters.” Journal of Mathematical Psychology 52 (1): 1–9.

  1. I think there may have been a sign error in the equation for \(g_j^+\) on the top of the right column of page 9: I think the argument inside \(\Phi\) should be negated. The corresponding argument in the Mills ratio further down the page appears to be correct.↩︎

  2. this can happen for the parameter values: rt = 30, response = 1, a = 5, v = 0, t0 = 0, w = 0.5, sv = 1.5, err_tol = 1e-6, method = "NCDF"↩︎

  3. Note that as \(x \to \infty\), \(1 - \Phi(x) = \Phi(-x) \to 0\) and \(\phi(x) \to 0\); this results in \(M(x) = \frac{0}{0}\). Similarly for \(x \to -\infty\), \(1 - \Phi(x) = \Phi(-x) \to 1\) and \(\phi(x) \to 0\); this results in \(M(x) = \frac{1}{0}\).↩︎