Speed up (#119)

Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: 27856297+dependabot-preview[bot]@users.noreply.github.com <27856297+dependabot-preview[bot]@users.noreply.github.com>

Author: 3 people

.pre-commit-config.yaml

@@ -3,7 +3,7 @@
 # R specific hooks: https://github.com/lorenzwalthert/precommit
 repos:
   - repo: https://github.com/lorenzwalthert/precommit
-    rev: v0.4.2
+    rev: v0.4.3
     hooks:
       - id: style-files
         args: [--style_pkg=styler, --style_fun=tidyverse_style]
@@ -91,6 +91,6 @@ repos:
         files: '\.Rhistory|\.RData|\.Rds|\.rds$'
         # `exclude: <regex>` to allow committing specific files.
   - repo: https://github.com/igorshubovych/markdownlint-cli
-    rev: v0.40.0
+    rev: v0.41.0
     hooks:
       - id: markdownlint

NAMESPACE

@@ -51,6 +51,7 @@ importFrom(stats,dbeta)
 importFrom(stats,dbinom)
 importFrom(stats,integrate)
 importFrom(stats,optimize)
+importFrom(stats,pbeta)
 importFrom(stats,quantile)
 importFrom(stats,rbeta)
 importFrom(stats,rbinom)

R/dbetabinom.R

@@ -124,14 +124,20 @@ dbetaMix <- function(x, par, weights, log = FALSE) {
   assert_numeric(weights, lower = 0, upper = 1, finite = TRUE, any.missing = FALSE)
   assert_true(all.equal(sum(weights), 1))
   assert_true(identical(length(weights), nrow(par)))
-  ret <- sum(weights * dbeta(x, par[, 1], par[, 2]))
+  degree <- length(weights)
+
+  component_densities <- matrix(
+    dbeta(rep(x, each = degree), par[, 1], par[, 2]),
+    nrow = degree,
+    ncol = length(x)
+  )
+  ret <- as.numeric(weights %*% component_densities)
   if (log) {
     log(ret)
   } else {
     ret
   }
 }
-dbetaMix <- Vectorize(dbetaMix, vectorize.args = "x")
 
 
 #' Beta-Mixture CDF
@@ -155,15 +161,25 @@ dbetaMix <- Vectorize(dbetaMix, vectorize.args = "x")
 #' @note `q` can be a vector.
 #'
 #' @example examples/pbetaMix.R
+#' @importFrom stats pbeta
 #' @export
 pbetaMix <- function(q, par, weights, lower.tail = TRUE) {
-  assert_number(q, lower = 0, upper = 1, finite = TRUE)
+  assert_numeric(q, lower = 0, upper = 1, finite = TRUE)
   assert_numeric(weights, lower = 0, upper = 1, finite = TRUE)
   assert_matrix(par)
   assert_flag(lower.tail)
-  sum(weights * pbeta(q, par[, 1], par[, 2], lower.tail = lower.tail))
+  .pbetaMix(q = q, par = par, weights = weights, lower.tail = lower.tail)
+}
+
+.pbetaMix <- function(q, par, weights, lower.tail) {
+  degree <- length(weights)
+  component_p <- matrix(
+    pbeta(rep(q, each = degree), par[, 1], par[, 2], lower.tail = lower.tail),
+    nrow = degree,
+    ncol = length(q)
+  )
+  as.numeric(weights %*% component_p)
 }
-pbetaMix <- Vectorize(pbetaMix, vectorize.args = "q")
 
 
 #' Beta-Mixture Quantile Function
@@ -186,23 +202,33 @@ pbetaMix <- Vectorize(pbetaMix, vectorize.args = "q")
 #' @example examples/qbetaMix.R
 #' @export
 qbetaMix <- function(p, par, weights, lower.tail = TRUE) {
-  f <- function(pi) {
-    pbetaMix(q = pi, par = par, weights = weights, lower.tail = lower.tail) - p
-  }
-  # Note: we give the lower and upper function values here in order to avoid problems for
-  # p = 0 or p = 1.
-  unirootResult <- uniroot(
-    f,
-    lower = 0, upper = 1,
-    f.lower = -p, f.upper = 1 - p,
-    tol = sqrt(.Machine$double.eps) # Increase the precision over default `tol`.
-  )
-  if (unirootResult$iter < 0) {
-    NA
-  } else {
-    assert_number(unirootResult$root)
-    assert_true(all.equal(f(unirootResult$root), 0, tolerance = 0.001))
-    unirootResult$root
-  }
+  assert_numeric(p, lower = 0, upper = 1)
+  assert_numeric(weights, lower = 0, upper = 1, finite = TRUE)
+  assert_matrix(par)
+  assert_flag(lower.tail)
+
+  grid <- seq(0, 1, len = 31)
+  f_grid <- .pbetaMix(grid, par, weights, lower.tail = lower.tail)
+
+  sapply(p, function(p) {
+    # special cases
+    if (p == 0) {
+      return(0)
+    }
+    if (p == 1) {
+      return(1)
+    }
+
+    diff <- f_grid - p
+    pos <- diff > 0
+    grid_interval <- c(grid[!pos][which.max(diff[!pos])], grid[pos][which.min(diff[pos])])
+
+    uniroot(
+      f = function(q) .pbetaMix(q, par, weights, lower.tail = lower.tail) - p,
+      interval = grid_interval,
+      f.lower = -p,
+      f.upper = 1 - p,
+      tol = sqrt(.Machine$double.eps)
+    )$root
+  })
 }
-qbetaMix <- Vectorize(qbetaMix, vectorize.args = "p")

R/postprob.R

@@ -81,18 +81,25 @@ postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALS
     if (missing(weights)) {
       weights <- rep(1, nrow(parE))
     }
-    betamixPost <- h_getBetamixPost(
-      x = x,
+    betamixPost <- lapply(
+      x,
+      h_getBetamixPost,
       n = n,
       par = parE,
       weights = weights
     )
+  } else {
+    assert_list(betamixPost)
+    assert_names(names(betamixPost), identical.to = c("par", "weights"))
+    betamixPost <- list(betamixPost)
   }
-  assert_list(betamixPost)
-  assert_names(names(betamixPost), identical.to = c("par", "weights"))
-  ret <- with(
+
+  ret <- vapply(
     betamixPost,
-    pbetaMix(q = p, par = par, weights = weights, lower.tail = FALSE)
+    FUN = function(bmp) {
+      .pbetaMix(q = p, par = bmp$par, weights = bmp$weights, lower.tail = FALSE)
+    },
+    FUN.VALUE = numeric(length(p))
   )
 
   if (log.p) {
@@ -101,4 +108,3 @@ postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALS
     ret
   }
 }
-postprob <- Vectorize(postprob, vectorize.args = "x")

R/postprobDist.R

@@ -164,7 +164,7 @@ postprobDist <- function(x,
   if (missing(weightsS)) {
     weightsS <- rep(1, nrow(parS))
   }
-  assert_number(n, lower = x, finite = TRUE)
+  assert_number(n, lower = min(x), finite = TRUE)
   assert_numeric(x, lower = 0, upper = n, finite = TRUE)
   assert_number(nS, lower = 0, finite = TRUE)
   assert_number(xS, lower = 0, upper = nS, finite = TRUE)
@@ -174,10 +174,10 @@ postprobDist <- function(x,
   assert_numeric(weightsS, lower = 0, finite = TRUE)
   assert_numeric(parE, lower = 0, finite = TRUE)
   assert_numeric(parS, lower = 0, finite = TRUE)
-  activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)
+
+  activeBetamixPost <- lapply(x, function(x) h_getBetamixPost(x = x, n = n, par = parE, weights = weights))
+
   controlBetamixPost <- h_getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
-  assert_names(names(activeBetamixPost), identical.to = c("par", "weights"))
-  assert_names(names(controlBetamixPost), identical.to = c("par", "weights"))
   if (relativeDelta) {
     epsilon <- .Machine$double.xmin
     integrand <- h_integrand_relDelta
@@ -186,26 +186,26 @@ postprobDist <- function(x,
     integrand <- h_integrand
   }
   bounds <- h_get_bounds(controlBetamixPost = controlBetamixPost)
-  intRes <- integrate(
-    f = integrand,
-    lower =
-      max(
-        bounds[1],
-        ifelse(relativeDelta, 0, 0 - delta)
-      ),
-    upper =
-      min(
-        ifelse(relativeDelta, 1, 1 - delta),
-        bounds[2]
-      ),
-    delta = delta,
-    activeBetamixPost = activeBetamixPost,
-    controlBetamixPost = controlBetamixPost
+
+  integral_results <- lapply(
+    seq_along(x),
+    function(i, this_posterior = activeBetamixPost, input_x = x) {
+      intRes <- integrate(
+        f = integrand,
+        lower = max(bounds[1], ifelse(relativeDelta, 0, 0 - delta)),
+        upper = min(ifelse(relativeDelta, 1, 1 - delta), bounds[2]),
+        delta = delta,
+        activeBetamixPost = this_posterior[[i]],
+        controlBetamixPost = controlBetamixPost
+      )
+      if (intRes$message == "OK") {
+        intRes$value
+      } else {
+        warning("Integration failed for posterior based on x =", input_x[i], "\n", intRes$message)
+        NA_real_
+      }
+    }
   )
-  if (intRes$message == "OK") {
-    intRes$value
-  } else {
-    stop(intRes$message)
-  }
+
+  unlist(integral_results)
 }
-postprobDist <- Vectorize(postprobDist, vectorize.args = "x")

examples/postprob.R

@@ -18,3 +18,13 @@ postprob(
     ),
   weights = c(0.6, 0.4)
 )
+
+postprob(
+  x = 0:23, n = 23, p = 0.60,
+  par =
+    rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    ),
+  weights = c(0.6, 0.4)
+)

man/postprob.Rd

@@ -84,6 +84,23 @@ test_that("pbetaMix gives the correct number result with beta-mixture", {
   expect_equal(result, 0.4768404, tolerance = 1e-5)
 })
 
+test_that("pbetaMix works for edge cases", {
+  result_ushape <- pbetaMix(
+    q = c(0, 1),
+    par = rbind(c(0.2, 0.4), c(3, .3)),
+    weights = c(0.6, 0.4)
+  )
+  expect_equal(result_ushape, c(0, 1))
+
+  result_vshape <- pbetaMix(
+    q = c(0, 1),
+    par = rbind(c(9, 4), c(1, 1)),
+    weights = c(0.6, 0.4)
+  )
+  expect_equal(result_vshape, c(0, 1))
+})
+
+
 test_that("The complement of pbetaMix can be derived with a different lower.tail flag", {
   result <- pbetaMix(
     q = 0.3,
@@ -184,6 +201,32 @@ test_that("dbetaMix gives the correct result as dbeta", {
   expect_equal(result, result2, tolerance = 1e-4)
 })
 
+test_that("dbetaMix handles edge cases", {
+  result_inf <- dbetaMix(
+    x = c(0, 1), par = rbind(c(0.2, 0.4), c(1, 1)),
+    weights = c(0.6, 0.4)
+  )
+  expect_equal(result_inf, c(Inf, Inf))
+
+  result_finite <- dbetaMix(
+    x = c(0, 1), par = rbind(c(2, 4), c(1, 1)),
+    weights = c(0.6, 0.4)
+  )
+  expect_equal(result_finite, c(0.4, 0.4))
+
+  result_right <- dbetaMix(
+    x = c(0, 1), par = rbind(c(0, 4), c(1, 1)),
+    weights = c(0.6, 0.4)
+  )
+  expect_equal(result_right, c(Inf, 0.4))
+
+  result_right <- dbetaMix(
+    x = c(NA, 1), par = rbind(c(0, 4), c(1, 1)),
+    weights = c(0.6, 0.4)
+  )
+  expect_equal(result_right, c(NA, 0.4))
+})
+
 # h_getBetamixPost ----
 
 test_that("h_getBetamixPost gives the correct beta-mixture parameters", {

tests/testthat/test-dbetabinom.R

@@ -65,23 +65,28 @@ test_that("h_decision_two_predprobDist gives correct result and list", {
 
 test_that("ocPredprobDist gives correct result and list when relativeDelta = FALSE", {
   set.seed(1989)
-  result <- ocPredprobDist(
-    nnE = c(10, 20, 30),
-    truep = 0.40,
-    deltaE = 0.5,
-    deltaF = 0.5,
-    relativeDelta = FALSE,
-    tT = 0.6,
-    phiU = 0.80,
-    phiFu = 0.7,
-    parE = c(1, 1),
-    parS = c(5, 25),
-    weights = 1,
-    weightsS = 1,
-    sim = 50,
-    nnF = c(10, 20, 30),
-    wiggle = TRUE,
-    decision1 = FALSE
+  expect_warning(
+    {
+      result <- ocPredprobDist(
+        nnE = c(10, 20, 30),
+        truep = 0.40,
+        deltaE = 0.5,
+        deltaF = 0.5,
+        relativeDelta = FALSE,
+        tT = 0.6,
+        phiU = 0.80,
+        phiFu = 0.7,
+        parE = c(1, 1),
+        parS = c(5, 25),
+        weights = 1,
+        weightsS = 1,
+        sim = 50,
+        nnF = c(10, 20, 30),
+        wiggle = TRUE,
+        decision1 = FALSE
+      )
+    },
+    "achieve convergence"
   )
   result_sum <- sum(result$oc[5:7])
   expect_equal(result_sum, 1)

tests/testthat/test-ocPredprobDist.R

@@ -53,3 +53,18 @@ test_that("postprob gives incrementally higher values with increased x", {
   )
   expect_true(is_lower < is_higher)
 })
+
+test_that("postprob works with vector x", {
+  result <- postprob(x = 0:23, n = 23, p = 0.60, par = c(0.6, 0.4))
+  expected <- c(
+    1.12066620085448e-10, 6.73786529927603e-09, 1.45879637562279e-07,
+    1.86374536434781e-06, 1.64656040420248e-05, 0.000108838231763851,
+    0.000564103325535708, 0.00236446983272442, 0.00819197194809839,
+    0.0238449136766029, 0.0590640325657381, 0.125847456119664,
+    0.232931221473374, 0.378259188739121, 0.54495891589689,
+    0.705949748288983, 0.835980805221058, 0.922929283049132,
+    0.970355725500809, 0.991009176245894, 0.997963909660055,
+    0.999685712592687, 0.999972679748126, 0.99999934483779
+  )
+  expect_equal(result, expected, tolerance = 1e-5)
+})