Benchmarking the Density Function Approximations
We want to determine the performance of each implementation across a variety of parameters in order to identify any slow areas for individual implementations. To achieve this rigor, we define a parameter space (in a code chunk below) and loop through each combination of parameters. For the preliminary benchmark testing, we will input the response times as a vector to each function instead of inputting individual response times for two reasons: first, this is the most common way to input the response time data; and second, this allows for the R-based implementation to exploit vectorization to make that implementation as fast as possible. For each combination of parameters, 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.
Upon collecting the benchmark data, we perform a preliminary analysis of the implementations’ performance by showing the distributions of benchmark times as side-by-side violin plots for each implementation. Then we cull the implementations and only test further on the fastest and most widely used implementations. To elucidate the areas of the parameter space where the different implementations excel and struggle, we plot the median benchmark times as a function of the response time input. This plot will reveal the efficiencies of the various implementations that we are testing.
Generating Benchmark Data
Before analyzing how the different density function approximations perform, we must first generate the benchmark data. In addition to testing all of the implementations available in dfddm(), we also include three functions that are considered to be current and widely used. First, we include the function rtdists::ddiffusion() from the package rtdists as it is well known for being not only user friendly and feature-rich but also designed specifically for handling data with regard to distributions for response time models. Second, we also test the RWiener::dwiener() function from the RWiener package, which is mainly aimed at providing an R language interface to calculate various functions from the Wiener process. Third, we include some raw R code that is bundled with the paper written by Gondan, Blurton, and Kesselmeier (2014) about improving the approximation to the density function. As this code was not yet available in an R package, it is part of fddm and can be accessed after running source on the corresponding file as shown below.
To capture the behavior of the density function approximations, we test each one in a granular parameter space that includes the parameter values most likely used in practice. This parameter space can be found fully defined in the code chunk later in this section. Note that we are only testing the “lower” density function as the conversion from the “lower” density function to the “upper” density function involves only the mappings \(v \to -v\) and \(w \to 1-w\). Since our parameter space includes the negatives of all positive values of \(v\) and the complements of all values of \(w\), testing both the “lower” and “upper” density function approximations would be redundant.
In addition, we include two functions for benchmarking the performance of the implementations: one where the response times are input as a vector, and one where the response times are input one at a time to the density function. Inputting the response times as a vector is the most common way of inputting data to the density function and allows the R code from Gondan, Blurton, and Kesselmeier (2014) to exploit R’s powerful vectorization, operating at its highest efficiency. We also include a benchmark function that inputs the response times individually as this allows us to see the impact that the size of the input response time has on the evaluation time of the density function approximations.
At each point in this parameter space, we run the microbenchmark::microbenchmark() function 1000 times for each implementation and save the median of these 1000 benchmark times. Using the median of these data should reduce the number of outliers in the data by limiting the influence of computational noise on the results. Note that for expositional purposes we will only say “benchmark times” to mean the collection of medians of the 1000 microbenchmark::microbenchmark() measurements recorded at each point in the parameter space. The following code chunk defines the functions that we will use to benchmark the implementations; we will be using the microbenchmark package to time the evaluations.
library("fddm")
library("rtdists")
library("RWiener")
source(system.file("extdata", "Gondan_et_al_density.R",
package = "fddm", mustWork = TRUE))
library("microbenchmark")
rt_benchmark_vec <- function(RT, resp, V, A, t0 = 1e-4, W = 0.5,
err_tol = 1e-6, times = 100, unit = "ns") {
fnames <- c("fs_SWSE_17", "fs_SWSE_14", "ft_SWSE_17", "ft_SWSE_14",
"fb_SWSE_17", "fb_SWSE_14",
"fs_Gon_17", "fs_Gon_14", "fb_Gon_17", "fb_Gon_14",
"fs_Nav_17", "fs_Nav_14", "fb_Nav_17", "fb_Nav_14",
"fl_Nav_09", "RWiener", "Gondan", "rtdists")
nf <- length(fnames) # number of functions being benchmarked
nV <- length(V)
nA <- length(A)
nW <- length(W)
resp <- rep(resp, length(RT)) # for RWiener
# Initialize the dataframe to contain the microbenchmark results
mbm_res <- data.frame(matrix(ncol = 3+nf, nrow = nV*nA*nW))
colnames(mbm_res) <- c('V', 'A', 'W', 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) {
mbm <- microbenchmark(
fs_SWSE_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "SWSE",
summation_small = "2017"),
fs_SWSE_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "SWSE",
summation_small = "2014"),
ft_SWSE_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "eff_rt", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2017"),
ft_SWSE_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "eff_rt", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2014"),
fb_SWSE_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms_large", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2017"),
fb_SWSE_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms_large", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2014"),
fs_Gon_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Gondan",
summation_small = "2017"),
fs_Gon_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Gondan",
summation_small = "2014"),
fb_Gon_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Gondan",
summation_small = "2017"),
fb_Gon_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Gondan",
summation_small = "2014"),
fs_Nav_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Navarro",
summation_small = "2017"),
fs_Nav_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Navarro",
summation_small = "2014"),
fb_Nav_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Navarro",
summation_small = "2017"),
fb_Nav_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Navarro",
summation_small = "2014"),
fl_Nav_09 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "large"),
RWiener = dwiener(RT, resp = resp, alpha = A[a],
delta = V[v], tau = t0, beta = W[w],
give_log = FALSE),
Gondan = fs(t = RT-t0, a = A[a], v = V[v],
w = W[w], eps = err_tol), # only "lower" resp
rtdists = ddiffusion(RT, resp, a = A[a], v = V[v],
t0 = t0, z = W[w]*A[a]),
times = times, unit = unit)
# 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]
# add the median microbenchmark results to the dataframe
for (i in 1:nf) {
mbm_res[row_idx, 3+i] <- median(mbm[mbm[,1] == fnames[i],2])
}
# iterate start value
row_idx <- row_idx + 1
}
}
}
return(mbm_res)
}
rt_benchmark_ind <- function(RT, resp, V, A, t0 = 1e-4, W = 0.5,
err_tol = 1e-6, times = 100, unit = "ns") {
fnames <- c("fs_SWSE_17", "fs_SWSE_14", "ft_SWSE_17", "ft_SWSE_14",
"fb_SWSE_17", "fb_SWSE_14",
"fs_Gon_17", "fs_Gon_14", "fb_Gon_17", "fb_Gon_14",
"fs_Nav_17", "fs_Nav_14", "fb_Nav_17", "fb_Nav_14",
"fl_Nav_09", "RWiener", "Gondan", "rtdists")
nf <- length(fnames) # number of functions being benchmarked
nRT <- length(RT)
nV <- length(V)
nA <- length(A)
nW <- length(W)
# Initialize the dataframe to contain the microbenchmark results
mbm_res <- data.frame(matrix(ncol = 4+nf, nrow = nRT*nV*nA*nW))
colnames(mbm_res) <- c('RT', 'V', 'A', 'W', 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) {
mbm <- microbenchmark(
fs_SWSE_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "SWSE",
summation_small = "2017"),
fs_SWSE_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "SWSE",
summation_small = "2014"),
ft_SWSE_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "eff_rt", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2017"),
ft_SWSE_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "eff_rt", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2014"),
fb_SWSE_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms_large", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2017"),
fb_SWSE_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms_large", switch_thresh = 0.8,
n_terms_small = "SWSE", summation_small = "2014"),
fs_Gon_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Gondan",
summation_small = "2017"),
fs_Gon_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Gondan",
summation_small = "2014"),
fb_Gon_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Gondan",
summation_small = "2017"),
fb_Gon_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Gondan",
summation_small = "2014"),
fs_Nav_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Navarro",
summation_small = "2017"),
fs_Nav_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "small", n_terms_small = "Navarro",
summation_small = "2014"),
fb_Nav_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Navarro",
summation_small = "2017"),
fb_Nav_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "terms", n_terms_small = "Navarro",
summation_small = "2014"),
fl_Nav_09 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
switch_mech = "large"),
RWiener = dwiener(RT[rt], resp = resp, alpha = A[a],
delta = V[v], tau = t0, beta = W[w],
give_log = FALSE),
Gondan = fs(t = RT[rt]-t0, a = A[a], v = V[v],
w = W[w], eps = err_tol), # only "lower" resp
rtdists = ddiffusion(RT[rt], resp, a = A[a], v = V[v],
t0 = t0, z = W[w]*A[a]),
times = times, unit = unit)
# 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]
# 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)
}Control loops in R have a lot of overhead and are comparatively slow compared to loops in C++; thus, we will not run rt_benchmark_vec() in this vignette, and we will load pre-run benchmark data instead. The pre-run data was collected using the same functions and parameter space as defined in this vignette. Since this benchmark function contains one fewer loop because of the vectorized response times, we increase the number of samples that the microbenchmark::microbenchmark() function considers when timing the various implementations. To show comparisons across the implementations with respect to the input response time, we will again load pre-run benchmark data obtained through using rt_benchmark_ind() as this takes even longer to evaluate than the benchmark with vectorized response times. The following code chunk defines the parameter space to be explored throughout the benchmark testing and shows an example of how to run the benchmark functions defined in the previous code chunk.
# 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 <- 1e-4 # must be nonzero for RWiener
W <- seq(0.3, 0.7, by = 0.1)
err_tol <- 1e-6 # this is the setting from rtdists
# Run benchmark tests
bm_vec <- rt_benchmark_vec(RT = RT, resp = "lower", V = V, A = A, t0 = t0,
W = W, err_tol = err_tol,
times = 1000, unit = "ns")
bm_ind <- rt_benchmark_ind(RT = RT, resp = "lower", V = V, A = A, t0 = t0,
W = W, err_tol = err_tol,
times = 100, unit = "ns")Analysis of Benchmark Results
The first step in analyzing the benchmark data is a rough visualization of the results in side-by-side violin plots that show the distribution of microbenchmark::microbenchmark() timings. We use the data that we just generated in the previous code chunk to plot a density violin for each implementation. Each density violin is overlaid with a boxplot showing the median and the first and third quartiles, in addition to a dashed line representing the mean. We also zoom in to the area of the plot where most of the fddm implementations lie; depending on your system, this may obscure the results from rtdists as they are typically much higher than that of fddm.
For this visualization we will use the benchmark data where the response times were handed over as a vector to each fddm call. As discussed above, running the benchmark on individual response times is considerably slower as the loop over the response times resides in R and not in C++.
library("reshape2")
library("ggplot2")
library("ggforce")
# load data, will be in the variable 'bm_vec'
load(system.file("extdata", "dfddm_density", "bm_vec_0-2.Rds",
package = "fddm", mustWork = TRUE))
t_idx <- match("W", 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("ft_SWSE_17", "ft_SWSE_14", "fb_SWSE_17", "fb_SWSE_14",
"fb_Gon_17", "fb_Gon_14", "fb_Nav_17", "fb_Nav_14",
"fs_SWSE_17", "fs_SWSE_14",
"fs_Gon_17", "fs_Gon_14", "fs_Nav_17", "fs_Nav_14",
"fl_Nav_09", "RWiener", "Gondan", "rtdists")
Color_vec <- c("#92c639", "#d3e8b0", "#729438", "#d2e3b5",
"#b3724d", "#e0c7b8", "#4da7b3", "#b8dce0",
"#5cc639", "#bee8b0",
"#b34d4d", "#e0b8b8", "#4d80b3", "#b8cce0",
"#dcdca3", "#deccba", "#c5a687", "#ac8053")
Outline_vec <- c("#92c639", "#92c639", "#729438", "#729438",
"#b3724d", "#b3724d", "#4da7b3", "#4da7b3",
"#5cc639", "#5cc639",
"#b34d4d", "#b34d4d", "#4d80b3", "#4d80b3",
"#dcdca3", "#deccba", "#c5a687", "#ac8053")
mi <- min(bm_vec[, -seq_len(t_idx)])
ma <- max(bm_vec[, (t_idx+1):(ncol(bm_vec)-4)])
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 = "none") +
scale_fill_manual(values = Color_vec, guide = "none") +
geom_boxplot(width = 0.15, fill = "white", alpha = 0.5) +
stat_summary(fun = mean, geom = "errorbar",
aes(ymax = after_stat(y), ymin = after_stat(y)),
width = .35, linetype = "dashed") +
scale_x_discrete(labels = c(
bquote(f[t] ~ SWSE[17]), bquote(f[t] ~ SWSE[14]),
bquote(f[c] ~ SWSE[17]), bquote(f[c] ~ SWSE[14]),
bquote(f[c] ~ Gon[17]), bquote(f[c] ~ Gon[14]),
bquote(f[c] ~ Nav[17]), bquote(f[c] ~ Nav[14]),
bquote(f[s] ~ SWSE[17]), bquote(f[s] ~ SWSE[14]),
bquote(f[s] ~ Gon[17]), bquote(f[s] ~ Gon[14]),
bquote(f[s] ~ Nav[17]), bquote(f[s] ~ Nav[14]),
bquote(f[l] ~ "Nav"), "RWiener", "Gondan", "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 small-time implementations (\(f_s\)) have long tails due to very few outliers, but their central tendencies are rather fast. The implementations combining small-time and large-time (\(f_c\)) approximations show an almost normal distribution with only a short tail; they also have rather fast central tendencies. Of the fddm implementations, the large-time implementation (\(f_\ell\)) is the slowest with a long tail. The function from the RWiener package combines small-time and large-time approach and also shows an almost normal distribution, but its distribution of median benchmark times is shifted higher compared to those of fddm. The other two implementations are clearly slower; the function provided by Gondan, Blurton, and Kesselmeier (2014) is purely in R so it is constrained by the speed of R’s vectorization. Similarly constrained by R, the function from the rtdists package performs quite a few computations in R that make it very user-friendly but at the cost of having the slowest central tendencies.
To visualize how the response time input can affect the efficiency of each implementation, we plot the benchmark times as a function of the effective response time. We define the effective response time as \(\frac{t}{a^2}\) because the two model parameters \(t\) and \(a\) always appear inside the infinite sums in this relationship. In the following plot, we will vary the effective response time and show the mean of the benchmark times, the 10% quantiles and 90% quantiles of the benchmark times, as well as the minimum and maximum benchmark times. For each value of the effective response time, we aggregate the benchmark data for all values of the other parameter to generate the statistics plotted. Note that because we vary the response time in this plot, we use the benchmark data generated from rt_benchmark_ind). As this benchmark function records the evaluation time of approximating the density for only one response time in contrast to a vector of response times, we expect the typical benchmark times produced by rt_benchmark_ind to be faster than that of rt_benchmark_vec.
To avoid too many repetitious plots, we will only include six of the implementations from the previous plot that illustrate the different variations and implementations of the density function approximations. We include: the fastest small-time implementation in dfddm() (“\(f_s \text{SWSE}_{17}\)”), the one large-time implementation of the large-time in dfddm (“\(f_\ell\) Nav”), the fastest implementation in dfddm() that combines the small-time and large-time (“\(f_c \text{SWSE}_{17}\)”), and the three implementations from R packages and the literature as described above.
The following code chunk loads the pre-run benchmark data, prepares it, and plots the benchmark times as the effective response time varies.
# load data, will be in the variable 'bm_ind'
load(system.file("extdata", "dfddm_density", "bm_ind_0-2.Rds",
package = "fddm", mustWork = TRUE))
bm_ind[["RTAA"]] <- bm_ind[["RT"]] / bm_ind[["A"]] / bm_ind[["A"]]
bm_ind <- bm_ind[, c(1, 2, ncol(bm_ind), 3:(ncol(bm_ind)-1)) ]
t_idx <- match("W", 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("ft_SWSE_17", "fs_SWSE_17", "fl_Nav_09",
"RWiener", "Gondan", "rtdists")
Color_meq <- c("#92c639", "#5cc639", "#dcdca3",
"#deccba", "#c5a687", "#ac8053")
mbm_meq <- subset(mbm_ind, FuncName %in% Names_meq)
my_labeller <- as_labeller(c(ft_SWSE_17 = "f[t] ~ SWSE[17]",
fs_SWSE_17 = "f[s] ~ SWSE[17]",
fl_Nav_09 = "f[l] ~ Nav",
RWiener = "RWiener",
Gondan = "Gondan",
rtdists = "rtdists"),
default = label_parsed)
ggplot(mbm_meq, aes(x = RTAA, 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.001, 0.1, 10, 1000),
labels = as.character(c(0.001, 0.1, 10, 1000))) +
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(frac(t, a^2) ~ ", effective response time, " ~ log[10]),
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 plot confirms that the aptly named small-time and large-time density functions are generally more efficient at calculating the density for small response times and large response times, respectively. Of the three implementations used from the literature, the RWiener and rtdists functions use a switching mechanism between the small-time and large-time density functions, and the implementation from Gondan, Blurton, and Kesselmeier (2014) only uses the small-time density function. Importantly, the implementations that utilize a switching mechanism between the small-time and large-time density function approximations have much flatter plots and generally perform well across the \(\frac{t}{a^2}\) parameter space; this confirms that using a combination of both timescales results in a very stable implementation.
