Speed up

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 @@
   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
@@ -150,20 +156,34 @@
 #' @typed lower.tail : flag
 #'  if `TRUE` (default), probabilities are `P[X <= x]`,
 #'  and otherwise `P[X > x]`.
+#' @typed skipchecks : flag
+#'   Don't check arguments `q`, `weights`, `par`, `lower.tail`.
+#'   Useful to speed up execution within functions where values
+#'   are known to be correct.
 #' @return The (one minus) cdf value
 #'
 #' @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(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 <- function(q, par, weights, lower.tail = TRUE, skipchecks = FALSE) {
+  if (isFALSE(skipchecks)) {
+    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)
+  }
+
+  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 +206,34 @@
 #' @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, skipchecks = TRUE)
+
+  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, skipchecks = TRUE) - p,
+      # interval = c(0, 1),
+      interval = grid_interval,
+      f.lower = -p,
+      f.upper = 1 - p,
+      tol = sqrt(.Machine$double.eps)
+    )$root
+  })
 }
-qbetaMix <- Vectorize(qbetaMix, vectorize.args = "p")

R/plotBeta.R

@@ -2,7 +2,6 @@
 #'
 #' This function will plot the PDF of a beta distribution
 #'
-#' @inheritParams dbetabinom
 #' @typed alpha : number
 #'  first parameter of the Beta distribution
 #' @typed beta : number

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, skipchecks = TRUE)
+    },
+    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 @@
   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,12 @@
   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"))
+  # 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 +188,26 @@
     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)
+)

tests/testthat/test-ocPredprobDist.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)