tests/testthat/test-partial-derivatives.R

In fddm: Fast Implementation of the Diffusion Decision Model

context("Testing the partial derivatives of the density function for accuracy")
# note: allowance is 4*err because WienR doesn't account for 2 sums

test_that("Accuracy relative to WienR", {
  testthat::skip_if_not_installed("WienR")
  library("WienR")

  ### Define Parameter Space ###
  RT <- c(0.001, 0.1, 1, 2, 3, 4, 5, 10, 30)
  A <- c(0.25, 0.5, 1, 2.5, 5)
  V <- c(-5, -2, 0, 2, 5)
  W <- c(0.2, 0.5, 0.8)
  SV <- c(0, 0.5, 1, 1.5)
  t0 <- 0
  err <- 1e-12 # default setting from WienR

  nRT <- length(RT)
  nA <- length(A)
  nV <- length(V)
  nW <- length(W)
  nSV <- length(SV)
  N <- nRT * nA * nV * nW * nSV

  resp <- "lower"
  respn <- rep(resp, N)
  rt <- rep(RT, each = nSV * nW * nV * nA, times = 1)
  a <-  rep(A,  each = nSV * nW * nV, times = nRT)
  v <-  rep(V,  each = nSV * nW, times = nRT * nA)
  w <-  rep(W,  each = nSV, times = nRT * nA * nV)
  sv <- rep(SV, each = 1, times = nRT * nA * nV * nW)



  ###
  ### dv ###
  df <- data.frame(
    rt = rt - t0,
    a = a,
    v = v,
    w = w,
    sv = sv,
    res_WienR = dvWienerPDF(t = rt, response = respn, a = a, v = v, w = w,
                            t0 = t0, sv = sv)[["deriv"]],
    res_fddm = dv_dfddm(rt = rt, response = resp, a = a, v = v, t0 = t0,
                        w = w, sv = sv, err_tol = err),
    diff = numeric(N)
  )
  df[["diff"]] <- df[["res_WienR"]] - df[["res_fddm"]]

  # agreement with WienR
  expect_true(all(abs(df[["diff"]]) <= 4*err))



  ###
  ### da ###
  df[["res_WienR"]] <- daWienerPDF(t = rt, response = respn, a = a, v = v,
                                   w = w, t0 = t0, sv = sv)[["deriv"]]
  df[["res_fddm"]] <- da_dfddm(rt = rt, response = resp, a = a, v = v, t0 = t0,
                               w = w, sv = sv, err_tol = err)
  df[["diff"]] <- df[["res_WienR"]] - df[["res_fddm"]]

  # agreement with WienR
  expect_true(sum(abs(df[["diff"]]) <= 4*err) / nrow(df) > 0.98)

  # disagreements with WienR
  df_bad <- df[abs(df[["diff"]]) > 4*err, ]

  # small-time and large-time in fddm are equal
  fddm_sm <- da_dfddm(rt = df_bad[["rt"]], response = resp,
                      a = df_bad[["a"]], v = df_bad[["v"]], w = df_bad[["w"]],
                      sv = df_bad[["sv"]], t0 = 0, sl_thresh = Inf,
                      err_tol = err)
  fddm_lg <- da_dfddm(rt = df_bad[["rt"]], response = resp,
                      a = df_bad[["a"]], v = df_bad[["v"]], w = df_bad[["w"]],
                      sv = df_bad[["sv"]], t0 = 0, sl_thresh = -1,
                      err_tol = err)
  expect_true(sum(abs(fddm_sm - fddm_lg) <= 2*err) / nrow(df_bad) == 1)

  # disagreements are small relative to the calculation
  # WienR precision not guaranteed for sv > 0
  df_bad <- df_bad[df_bad[["sv"]] == 0, ]
  expect_true(all(abs(df_bad[["diff"]]) / abs(df_bad[["res_WienR"]]) <=
                  err*4e-2))



  ###
  ### dt ###
  df[["res_WienR"]] <- dtWienerPDF(t = rt, response = respn, a = a, v = v,
                                   w = w, t0 = t0, sv = sv)[["deriv"]]
  df[["res_fddm"]] <- dt_dfddm(rt = rt, response = resp, a = a, v = v, t0 = t0,
                               w = w, sv = sv, err_tol = err)
  df[["diff"]] <- df[["res_WienR"]] - df[["res_fddm"]]

  # agreement with WienR
  expect_true(sum(abs(df[["diff"]]) <= 4*err) / nrow(df) > 0.97)

  # disagreements with WienR (WienR precision not guaranteed for sv > 0)
  df_bad <- df[abs(df[["diff"]]) > 4*err & df[["sv"]] == 0, ]

  # disagreements are small relative to the calculation
  expect_true(all(abs(df_bad[["diff"]]) / abs(df_bad[["res_WienR"]]) <=
                  err*4e-2))



  ###
  ### dt0 ###
  df[["res_WienR"]] <- dt0WienerPDF(t = rt, response = respn, a = a, v = v,
                                    w = w, t0 = t0, sv = sv)[["deriv"]]
  df[["res_fddm"]] <- dt0_dfddm(rt = rt, response = resp, a = a, v = v, t0 = t0,
                                w = w, sv = sv, err_tol = err)
  df[["diff"]] <- df[["res_WienR"]] - df[["res_fddm"]]

  # agreement with WienR
  expect_true(sum(abs(df[["diff"]]) <= 4*err) / nrow(df) > 0.97)

  # disagreements with WienR (WienR precision not guaranteed for sv > 0)
  df_bad <- df[abs(df[["diff"]]) > 4*err & df[["sv"]] == 0, ]

  # disagreements are small relative to the calculation
  expect_true(all(abs(df_bad[["diff"]]) / abs(df_bad[["res_WienR"]]) <=
                  err*4e-2))



  ###
  ### dw ###
  df[["res_WienR"]] <- dwWienerPDF(t = rt, response = respn, a = a, v = v,
                                   w = w, t0 = t0, sv = sv)[["deriv"]]
  df[["res_fddm"]] <- dw_dfddm(rt = rt, response = resp, a = a, v = v, t0 = t0,
                               w = w, sv = sv, err_tol = err)
  df[["diff"]] <- df[["res_WienR"]] - df[["res_fddm"]]

  # agreement with WienR
  expect_true(sum(abs(df[["diff"]]) <= 4*err) / nrow(df) > 0.98)

  # disagreements with WienR
  df_bad <- df[abs(df[["diff"]]) > 4*err, ]

  # small-time and large-time in fddm are equal
  fddm_sm <- dw_dfddm(rt = df_bad[["rt"]], response = resp,
                      a = df_bad[["a"]], v = df_bad[["v"]], w = df_bad[["w"]],
                      sv = df_bad[["sv"]], t0 = 0, sl_thresh = -1,
                      err_tol = err)
  fddm_lg <- dw_dfddm(rt = df_bad[["rt"]], response = resp,
                      a = df_bad[["a"]], v = df_bad[["v"]], w = df_bad[["w"]],
                      sv = df_bad[["sv"]], t0 = 0, sl_thresh = Inf,
                      err_tol = err)
  expect_true(sum(abs(fddm_sm - fddm_lg) <= 2*err) / nrow(df_bad) == 1)



  ###
  ### dsv ###
  SV <- c(0.5, 1, 1.5) # remove sv = 0
  nSV <- length(SV)
  N <- nRT * nA * nV * nW * nSV
  rt <- rep(RT, each = nSV * nW * nV * nA, times = 1) - t0
  a <-  rep(A,  each = nSV * nW * nV, times = nRT)
  v <-  rep(V,  each = nSV * nW, times = nRT * nA)
  w <-  rep(W,  each = nSV, times = nRT * nA * nV)
  sv <- rep(SV, each = 1, times = nRT * nA * nV * nW)
  N <- nRT * nA * nV * nW * nSV
  respn <- rep(resp, N)

  df <- data.frame(
    rt = rt - t0,
    a = a,
    v = v,
    w = w,
    sv = sv,
    res_WienR = dsvWienerPDF(t = rt, response = respn, a = a, v = v,
                             w = w, t0 = t0, sv = sv)[["deriv"]],
    res_fddm = dsv_dfddm(rt = rt, response = resp, a = a, v = v, t0 = t0,
                         w = w, sv = sv, err_tol = err),
    diff = numeric(N)
  )
  df[["diff"]] <- df[["res_WienR"]] - df[["res_fddm"]]

  # agreement with WienR
  expect_true(all(abs(df[["diff"]]) <= 4*err))
})





test_that("Accuracy relative to numerical approximations", {
  testthat::skip_if_not_installed("numDeriv")
  library("numDeriv")

  ### Define Parameter Space ###
  RT <- c(0.001, 0.1, 1, 2, 3, 4, 5, 10, 30)
  A <- c(0.25, 0.5, 1, 2.5, 5)
  V <- c(-5, -2, 0, 2, 5)
  W <- c(0.2, 0.5, 0.8)
  SV <- c(0, 0.5, 1, 1.5)
  t0 <- 0
  err <- 1e-12
  eps <- 2 * sqrt(err)

  nRT <- length(RT)
  nA <- length(A)
  nV <- length(V)
  nW <- length(W)
  nSV <- length(SV)
  N <- nRT * nA * nV * nW * nSV

  resp <- "lower"
  rt <- rep(RT, each = nSV * nW * nV * nA, times = 1) - t0
  a <-  rep(A,  each = nSV * nW * nV, times = nRT)
  v <-  rep(V,  each = nSV * nW, times = nRT * nA)
  w <-  rep(W,  each = nSV, times = nRT * nA * nV)
  sv <- rep(SV, each = 1, times = nRT * nA * nV * nW)



  ###
  ### dv ###
  df <- data.frame(
    rt = rt,
    a = a,
    rtaa = rt/(a*a),
    v = v,
    w = w,
    sv = sv,
    res_small = dv_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                         w = w, sv = sv, sl_thresh = 1000, err_tol = err),
    res_large = dv_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                         w = w, sv = sv, sl_thresh = 0, err_tol = err),
    res_appx = numeric(N),
    diff_small = numeric(N),
    diff_large = numeric(N)
  )
  dv_wrap <- function(v, p) {
    return(dfddm(rt = p[1], response = resp, v = v, a = p[2], t0 = 0,
                 w = p[3], sv = p[4], err_tol = err))
  }
  for (i in seq_len(N)) {
    df[i, "res_appx"] <-
      numDeriv::grad(func = dv_wrap, x = c("v" = v[i]), method = "Richardson",
                    p = c(rt = rt[i], a = a[i], w = w[i], sv = sv[i]))
  }
  df[["diff_small"]] <- df[["res_appx"]] - df[["res_small"]]
  df[["diff_large"]] <- df[["res_appx"]] - df[["res_large"]]

  expect_true(all(abs(df[["diff_small"]]) <= eps))
  expect_true(all(abs(df[["diff_large"]]) <= eps))



  ###
  ### da ###
  df[["res_small"]] <- da_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                                w = w, sv = sv, sl_thresh = 1000, err_tol = err)
  df[["res_large"]] <- da_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                                w = w, sv = sv, sl_thresh = 0, err_tol = err)
  da_wrap <- function(a, p) {
    return(dfddm(rt = p[1], response = resp, v = p[2], a = a, t0 = 0,
                 w = p[3], sv = p[4], err_tol = err))
  }
  for (i in seq_len(N)) {
    df[i, "res_appx"] <-
      numDeriv::grad(func = da_wrap, x = c("a" = a[i]), method = "Richardson",
                     p = c(rt = rt[i], v = v[i], w = w[i], sv = sv[i]))
  }
  df[["diff_small"]] <- df[["res_appx"]] - df[["res_small"]]
  df[["diff_large"]] <- df[["res_appx"]] - df[["res_large"]]

  expect_true(all(abs(df[["diff_small"]]) <= eps))
  expect_true(all(abs(df[df[["rtaa"]] > 4e-03, "diff_large"]) <= eps))



  ###
  ### dt ###
  df[["res_small"]] <- dt_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                                w = w, sv = sv, sl_thresh = 1000, err_tol = err)
  df[["res_large"]] <- dt_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                                w = w, sv = sv, sl_thresh = 0, err_tol = err)
  dt_wrap <- function(rt, p) {
    return(dfddm(rt = rt, response = resp, v = p[1], a = p[2], t0 = 0,
                 w = p[3], sv = p[4], err_tol = err))
  }
  for (i in seq_len(N)) {
    df[i, "res_appx"] <-
      numDeriv::grad(func = dt_wrap, x = c("rt" = rt[i]), method = "Richardson",
                     p = c(v = v[i], a = a[i], w = w[i], sv = sv[i]))
  }
  df[["diff_small"]] <- df[["res_appx"]] - df[["res_small"]]
  df[["diff_large"]] <- df[["res_appx"]] - df[["res_large"]]

  expect_true(sum(abs(df[["diff_small"]]) <= eps) / N > 0.99)
  expect_true(sum(abs(df[["diff_large"]]) <= eps) / N > 0.98)

  # disagreements are small relative to the calculation
  expect_true(all(abs(df[abs(df[["diff_small"]]) > eps, "diff_small"]) /
                  abs(df[abs(df[["diff_small"]]) > eps, "res_appx"]) <= eps))
  expect_true(all(abs(df[abs(df[["diff_large"]]) > eps & df[["rtaa"]] > 0.01,
                         "diff_small"]) /
                  abs(df[abs(df[["diff_large"]]) > eps & df[["rtaa"]] > 0.01,
                         "res_appx"]) <= eps))



  ###
  ### dt0 ###
  df[["res_small"]] <- dt0_dfddm(rt = rt, response = resp, v = v, a = a,
                                 t0 = t0, w = w, sv = sv, sl_thresh = 1000,
                                 err_tol = err)
  df[["res_large"]] <- dt0_dfddm(rt = rt, response = resp, v = v, a = a,
                                 t0 = t0, w = w, sv = sv, sl_thresh = 0,
                                 err_tol = err)
  dt0_wrap <- function(rt, p) {
    return(dfddm(rt = rt, response = resp, v = p[1], a = p[2], t0 = 0,
                 w = p[3], sv = p[4], err_tol = err))
  }
  for (i in seq_len(N)) {
    df[i, "res_appx"] <- -1 *
      numDeriv::grad(func = dt0_wrap, x = c("rt" = rt[i]),
                     method = "Richardson",
                     p = c(v = v[i], a = a[i], w = w[i], sv = sv[i]))
  }
  df[["diff_small"]] <- df[["res_appx"]] - df[["res_small"]]
  df[["diff_large"]] <- df[["res_appx"]] - df[["res_large"]]

  expect_true(sum(abs(df[["diff_small"]]) <= eps) / N > 0.99)
  expect_true(sum(abs(df[["diff_large"]]) <= eps) / N > 0.98)

  # disagreements are small relative to the calculation
  expect_true(all(abs(df[abs(df[["diff_small"]]) > eps, "diff_small"]) /
                  abs(df[abs(df[["diff_small"]]) > eps, "res_appx"]) <= eps))
  expect_true(all(abs(df[abs(df[["diff_large"]]) > eps & df[["rtaa"]] > 0.01,
                         "diff_small"]) /
                  abs(df[abs(df[["diff_large"]]) > eps & df[["rtaa"]] > 0.01,
                         "res_appx"]) <= eps))



  ###
  ### dw ###
  df[["res_small"]] <- dw_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                                w = w, sv = sv, sl_thresh = 1000, err_tol = err)
  df[["res_large"]] <- dw_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                                w = w, sv = sv, sl_thresh = 0, err_tol = err)
  dw_wrap <- function(w, p) {
    return(dfddm(rt = p[1], response = resp, v = p[2], a = p[3], t0 = 0,
                 w = w, sv = p[4], err_tol = err))
  }
  for (i in seq_len(N)) {
    df[i, "res_appx"] <-
      numDeriv::grad(func = dw_wrap, x = c("w" = w[i]), method = "Richardson",
                    p = c(rt = rt[i], v = v[i], a = a[i], sv = sv[i]))
  }
  df[["diff_small"]] <- df[["res_appx"]] - df[["res_small"]]
  df[["diff_large"]] <- df[["res_appx"]] - df[["res_large"]]

  expect_true(all(abs(df[["diff_small"]]) <= eps))
  expect_true(all(abs(df[["diff_large"]]) <= eps))



  ###
  ### dsv2 ###
  SV <- c(0.5, 1, 1.5) # remove sv = 0
  nSV <- length(SV)
  N <- nRT * nA * nV * nW * nSV
  rt <- rep(RT, each = nSV * nW * nV * nA, times = 1) - t0
  a <-  rep(A,  each = nSV * nW * nV, times = nRT)
  v <-  rep(V,  each = nSV * nW, times = nRT * nA)
  w <-  rep(W,  each = nSV, times = nRT * nA * nV)
  sv <- rep(SV, each = 1, times = nRT * nA * nV * nW)

  dsv_wrap <- function(sv, p) {
    return(dfddm(rt = p[1], response = resp, v = p[2], a = p[3], t0 = 0,
                 w = p[4], sv = sv, err_tol = err))
  }
  df <- data.frame(
    rt = rt,
    a = a,
    v = v,
    w = w,
    sv = sv,
    res_small = dsv_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                          w = w, sv = sv, sl_thresh = 1000, err_tol = err),
    res_large = dsv_dfddm(rt = rt, response = resp, v = v, a = a, t0 = t0,
                          w = w, sv = sv, sl_thresh = 0, err_tol = err),
    res_appx = numeric(N),
    diff_small = numeric(N),
    diff_large = numeric(N)
  )
  for (i in seq_len(N)) {
    df[i, "res_appx"] <-
      numDeriv::grad(func = dsv_wrap, x = c("sv" = sv[i]),
                     method = "Richardson",
                     p = c(rt = rt[i], v = v[i], a = a[i], w = w[i]))
  }
  df[["diff_small"]] <- df[["res_appx"]] - df[["res_small"]]
  df[["diff_large"]] <- df[["res_appx"]] - df[["res_large"]]

  expect_true(all(abs(df[["diff_small"]]) <= eps))
  expect_true(all(abs(df[["diff_large"]]) <= eps))
})