Description of Methods in
pfddm
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 ∈ (0, ∞) (threshold separation),
v ∈ (−∞, ∞) (drift rate),
t0 ∈ [0, ∞)
(non-decision time/response time constant), w ∈ (0, 1) (relative starting
point), sv ∈ (0, ∞)
(inter-trial-variability of drift), and σ ∈ (0, ∞) (diffusion coefficient of
the underlying Wiener Process). Please note that for this vignette, we
will refer to the inter-trial variability of drift as η 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 = ∞. 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 Fupper(t | v, η, a, w) = Flower(t | − v, η, a, 1 − w).
Since the CDF of the DDM is widely used in parameter estimation
usually involving numerical optimization (e.g., χ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., η = 0) (see Blurton, Kesselmeier, and Gondan (2012) for details). There is no analytic large-time variant for η > 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 η > 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, ϵ, is bounded above by the absolute
value of the next term in the sum. For example, if the sum includes
k terms (b1 + b2 + … + bk),
then the approximation error is at most the absolute value of the k + 1th term (ϵ ≤ |bk + 1|).
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., η > 0), they first arrived at the following equation for the distribution function1:
where
and Φ(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 → ∞. 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
where Φ(x) is the standard normal CDF evaluated at x, and ϕ(x) is the standard normal PDF evaluated at x. Then the CDF of the DDM becomes
where again
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 (t0).
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 rtdistsNote 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 → −v and w → 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.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 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=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: Etc/UTC
#> 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.5.0 RWiener_1.3-3 rtdists_0.11-5
#> [7] fddm_1.0-2 rmarkdown_2.29
#>
#> 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.37 magrittr_2.0.3 estimability_1.5.1
#> [9] evaluate_1.0.1 grid_4.4.2 mvtnorm_1.3-2 fastmap_1.2.0
#> [13] plyr_1.8.9 jsonlite_1.8.9 Matrix_1.7-1 ggnewscale_0.5.0
#> [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.3
#> [25] crayon_1.5.3 rlang_1.1.4 expm_1.0-0 polyclip_1.10-7
#> [29] gsl_2.1-8 munsell_0.5.1 splines_4.4.2 withr_3.0.2
#> [33] cachem_1.1.0 yaml_2.3.10 tools_4.4.2 colorspace_2.1-1
#> [37] buildtools_1.0.0 vctrs_0.6.5 R6_2.5.1 emmeans_1.10.5
#> [41] lifecycle_1.0.4 stringr_1.5.1 MASS_7.3-61 pkgconfig_2.0.3
#> [45] pillar_1.9.0 bslib_0.8.0 gtable_0.3.6 glue_1.8.0
#> [49] Rcpp_1.0.13-1 tidyselect_1.2.1 xfun_0.49 tibble_3.2.1
#> [53] sys_3.4.3 knitr_1.49 msm_1.8.2 farver_2.1.2
#> [57] htmltools_0.5.8.1 labeling_0.4.3 maketools_1.3.1 compiler_4.4.2
#> [61] evd_2.3-7.1References
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/.
I think there may have been a sign error in the equation for gj+ on the top of the right column of page 9: I think the argument inside Φ should be negated. The corresponding argument in the Mills ratio further down the page appears to be correct.↩︎
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"↩︎Note that as x → ∞, 1 − Φ(x) = Φ(−x) → 0 and ϕ(x) → 0; this results in $M(x) = \frac{0}{0}$. Similarly for x → −∞, 1 − Φ(x) = Φ(−x) → 1 and ϕ(x) → 0; this results in $M(x) = \frac{1}{0}$.↩︎