Version 0.1.2

NAMESPACE

@@ -41,6 +41,7 @@ importFrom(scales,percent)
 importFrom(stats,density)
 importFrom(stats,dunif)
 importFrom(stats,dweibull)
+importFrom(stats,optimize)
 importFrom(stats,runif)
 importFrom(utils,capture.output)
 importFrom(utils,head)

NEWS.md

@@ -26,8 +26,6 @@
 
 # AcceptReject 0.1.2
 
-* The function `best_y()` has been implemented in Rcpp (C++) to enhance the candidate for the value of `y` that maximizes `f(y)/f_base(y)`. The `best_y()` function is internally called by the accept_reject() function and is not exported to the user. The `best_y()` function takes the form `best_y(NumericVector xlim, Function f, Function f_base, bool continuous, double epsilon = 0.001)` and returns a double;
-
 * The performance of the `one_step()` function, an internal function used in the implementation of C++ using Rcpp, has been improved.
 
 * The function `accept_reject()` now has the argument cores, which allows the user to control the number of cores that will be used if `parallel = TRUE`. The default, `cores = NULL`, means that all processor cores will be used. If `parallel = FALSE`, the cores argument is ignored.

R/RcppExports.R

@@ -5,7 +5,3 @@ one_step <- function(n, f, f_base, random_base, c) {
     .Call(`_AcceptReject_one_step`, n, f, f_base, random_base, c)
 }
 
-best_y <- function(xlim, f, f_base, continuous, epsilon = 0.001) {
-    .Call(`_AcceptReject_best_y`, xlim, f, f_base, continuous, epsilon)
-}
-

R/accept_reject.r

@@ -107,6 +107,7 @@
 #'
 #' @import rlang
 #' @importFrom lbfgs lbfgs
+#' @importFrom stats optimize
 #' @importFrom purrr partial map_dbl
 #' @importFrom numDeriv grad
 #' @importFrom parallel detectCores mclapply
@@ -169,33 +170,44 @@ accept_reject <-
     # Uniform distribution will be used if not all information from the base
     # distribution is provided.
     any_null <- any(is.null(c(f_base, random_base, args_f_base)))
+
     if (continuous && any_null) {
+      limit <- .Machine$double.xmin
+      if (xlim[1L] >= 0 && xlim[1L] <= limit) xlim[1L] <- limit
       f_base <- purrr::partial(.f = dunif, min = xlim[1L], max = xlim[2L])
       random_base <- purrr::partial(.f = runif, min = xlim[1L], max = xlim[2L])
     }
 
     # Is it a discrete random variable?
     if (!continuous) {
+      y <- xlim[1L]:xlim[2L]
       f_base <- function(x) dunif(x, min = xlim[1L], max = xlim[2L])
       random_base <- function(n) sample(x = xlim[1L]:xlim[2L], size = n, replace = TRUE)
-    } else if (continuous && !any_null) {
-      if (xlim[1L] == 0) xlim[1L] <- .Machine$double.xmin
+      ymax <- y[which.max(f(y)/f_base(y))]
+    } else {
       f_base <- purrr::partial(.f = f_base, !!!args_f_base)
       random_base <- purrr::partial(.f = random_base, !!!args_f_base)
-    }
 
-    ymax <- best_y(
-      xlim = xlim,
-      f = f,
-      f_base = f_base,
-      continuous = continuous,
-      epsilon = 0.001
-    )
+      objective_y <- function(y)
+        log(f_base(y)) - log(f(y))
 
-    objective_c <- function(c, y) {
+      try_objective_y <- function(...) {
+        tryCatch(
+          objective_y(...),
+          error = function(e) NaN
+        )
+      }
 
-      differences <- log(f(y)) - (log(c) + f_base(y))
+      suppressWarnings(ymax <-
+        optimize(
+          f = try_objective_y,
+          interval = xlim,
+        )$minimum
+      )
+    }
 
+    objective_c <- function(c, y) {
+      differences <- log(f(y)) - (log(c) + f_base(y))
       if (is.infinite(differences) && continuous) {
         return(.Machine$double.xmax)
       } else {
@@ -205,7 +217,7 @@ accept_reject <-
 
     gradient_objective_c <- function(c, y) {
       numDeriv::grad(
-        func = function(c) objective_c(c = c, y = ymax),
+        func = function(c) objective_c(c = c, y = y),
         x = c
       )
     }

README.md

@@ -284,7 +284,7 @@ case_1 <- accept_reject(
   xlim = c(0, 10)
 )
 toc()
-#> 0.385 sec elapsed
+#> 0.371 sec elapsed
 
 # Specifying the base probability density function
 tic()
@@ -300,7 +300,7 @@ case_2 <- accept_reject(
   c = 1.2
 )
 toc()
-#> 0.151 sec elapsed
+#> 0.145 sec elapsed
 
 # Visualizing the results
 p1 <- plot(case_1)