---
title: "Description of Methods in `pfddm`"
author: "Kendal Foster and Henrik Singmann"
date: '`r format(Sys.Date(), "%B %d, %Y")`'
bibliography: references.bib
output:
rmarkdown::html_vignette:
css: stile.css
toc: false
vignette: >
%\VignetteIndexEntry{Description of Methods in `pfddm`}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r echo=FALSE}
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)
}
```
```{r setup, include=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.
# Mathematical Background {#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](#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 {#dist}
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](math.html) 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 {#dist-inf}
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 @blurton2012fast for details). There is no analytic large-time variant for $\eta > 0$ because the resulting integral is incalculable (see @tuerlinckx2004efficient for more details). However, @blurton2017first provide the small-time variant of the CDF for $\eta > 0$; this variant contains an infinite sum that must be approximated. @blurton2017first 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](#bechmarks), 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 @blurton2017first.
### Normal CDF {#dist-inf-ncdf}
When @blurton2017first 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 function^[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.]:
\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 useful^[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"`].
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.
### Mills Ratio {#dist-inf-mills}
After deriving Equation $\ref{eq:cdf_ncdf}$, @blurton2017first manipulated their result to use the Mills ratio for its favorable numerical properties (see @blurton2012fast 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 ratio^[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}$.].
Fortunately, there exists a useful approximation to the Mills ratio for large arguments. @abramowitz1964handbook 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 [@sn] for `R`, which in turn was recommended by @gondan2014even. 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](#ben-res). 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 {#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`](https://cran.r-project.org/package=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](#par-space)) 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.
```{r bm-par-space}
# 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`](https://cran.r-project.org/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.
```{r bm-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 {#ben-vec}
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 @blurton2017first 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.
```{r bm-vec-fun}
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.
```{r bm-vec-run, eval=FALSE}
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)
```
```{r bm-vec-save, eval=FALSE, include=FALSE}
save(bm_vec, compress = "xz", compression_level = 9,
file = "../inst/extdata/pfddm_distribution/bm_vec_0-2.Rds")
```
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.
```{r bm-vec, fig.height=5, fig.width=8}
# 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")
```
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`](https://cran.r-project.org/package=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 @abramowitz1964handbook allows for a fast and accurate approximation to it.
## Non-vectorized Benchmark Data and Results {#ben-ind}
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 @blurton2017first 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.
```{r bm-ind-fun}
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.
```{r bm-ind-run, eval=FALSE}
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)
```
```{r bm-ind-save, eval=FALSE, include=FALSE}
save(bm_ind, compress = "xz", compression_level = 9,
file = "../inst/extdata/pfddm_distribution/bm_ind_0-2.Rds")
```
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.
```{r bm-meq-prep, fig.height=5, fig.width=8}
# 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](benchmark.html#den-ana). 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](#dist-inf), 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`](https://cran.r-project.org/package=rtdists); they are also faster than the implementation in the pure `R` code from @blurton2017first.
# {.unlisted .unnumbered}
#### R Session Info {.unlisted .unnumbered}
```{r session-info, collapse=TRUE}
sessionInfo()
```
# References