Error with rpowbend() - Polylogarithmic approximation

[AlbanSagouis]

Hi,

I'm working with mobiodiv/mobsim and I have errors with rpowbend (both with the CRAN version of sads and the latest version on the dev branch (bb004c1)):

> sads::rpowbend(n = 10, s = 2, omega = 0.5)
Error in if (cum + my.sq > want) break : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In (function (x, s, omega = 0.01, oM = -log(omega), log = FALSE)  :
  NaNs produced
> traceback()
7: qfinder(dpowbend, p[1], list(s = s, oM = oM), 0)
6: withCallingHandlers(expr, warning = function(w) if (inherits(w, 
       classes)) tryInvokeRestart("muffleWarning"))
5: suppressWarnings(qfinder(dpowbend, p[1], list(s = s, oM = oM), 
       0))
4: (function (p, s, omega = 0.01, oM = -log(omega), lower.tail = TRUE, 
       log.p = FALSE) 
   {
       if (!lower.tail) 
           p <- 1 - p
       op <- rank(p, ties.method = "max")
       p <- sort(p)
       if (length(s) > 1 | length(oM) > 1) 
           stop("Vectorization not implemented for the parameters")
       if (!missing(omega) && !missing(oM)) 
           stop("specify 'omega' or 'oM' but not both")
       if (log.p) 
           p <- exp(p)
       y <- c()
       y[1] <- suppressWarnings(qfinder(dpowbend, p[1], list(s = s, 
           oM = oM), 0))
       if (length(p) > 1) 
           for (i in 2:length(p)) y[i] <- suppressWarnings(qfinder(dpowbend, 
               p[i], list(s = s, oM = oM), y[i - 1]))
       if (any(is.nan(y))) 
    ...
3: do.call(qf, c(list(p = rr), coef))
2: shift_r("powbend", n, list(s = s, omega = omega))
1: sads::rpowbend(n = 10, s = 2, omega = 0.5)

Is it because of the argument values I use?

Thanks!

[piklprado]

Hi Alban

rpowbend relies on a approximation of the Polylogarithmic function that does not work for integer values of the parameter s larger than 1. To solve this limitation is to be done, but meanwhile you can circumvent the issue by specifying non-integer values for s >1. For instance:

> rpowbend(n = 10, s = 2.0001, omega = 0.5)
[1] 1 1 1 3 1 1 1 1 1 1

If you want to take a look at the function that approximates Polylog please see the function LiE at R/util.R . This function is not exported to the package workspace, as all other accessory functions.

[piklprado]

Thanks for poiting this isssue. I am keeping it open as a needed enhancement of the polylogarithmic approximation.