inst/doc/pfddm.R

In fddm: Fast Implementation of the Diffusion Decision Model

## ----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)
}

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  error = TRUE,
  comment = "#>"
)

## ----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

## ----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")

## ----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)
}

## ----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)

## ----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")

## ----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")

## ----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)
}

## ----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)

## ----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")

## ----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)

## ----session-info, collapse=TRUE----------------------------------------------
sessionInfo()