Why can't I use uniroot with plot?

Question

I am working on code that uses the uniroot function to approximate the root of an equation. I am trying to plot the behaviour of the function being passed through uniroot as the value of a free variable changes:

library(Deriv)

f1 <- function(s) {
  (1 - 2*s)^(-3/2)*exp((8*s)/(1-2*s))
}

f2 <- function(s) {
  log(f1(s))
}

f3 <- Deriv(f2, 's')
f4 <- Deriv(f3, 's')
f5 <- Deriv(f4, 's')

upp_s <- 1/2 - 1e-20

f_est <- function(x) {
  f3a <- function(s) {f3(s = s) - x}
  
  s_ <- uniroot(f3a,
                lower = -9,
                upper = upp_s)$root
  
  return(s_)
}

plot(f_est, from = 0, to=100, col="red", main="header")

The output of f_est works as expected. However, when passed through the plot function, uniroot seems to break:

> plot(f_est, from = 0, to=100, col="red", main="header")
Error in uniroot(f3a, lower = -9, upper = upp_s) : 
  f() values at end points not of opposite sign
In addition: Warning messages:
1: In if (is.na(f.lower)) stop("f.lower = f(lower) is NA") :
  the condition has length > 1 and only the first element will be used
2: In if (is.na(f.upper)) stop("f.upper = f(upper) is NA") :
 
 Error in uniroot(f3a, lower = -9, upper = upp_s) : 
f() values at end points not of opposite sign

The function is set up such that the endpoints specified in uniroot are always of opposite sign, and that there is always exactly one real root. I have also checked to confirm that the endpoints are non-missing when f_est is run by itself. I've tried vectorising the functions involved to no avail.

Why is this happening?

The code is giving `Error in Deriv(K_s, "s") : object 'K_s' not found`. — Rui Barradas, Sep 13 '22 at 16:07
I'm getting "f.upper = f(upper) is NA" for any value I put into `f_est()` (e.g. `f_est(1)`) — Ben Bolker, Sep 13 '22 at 17:26
@BenBolker This is because `upp_s` iss nowhere defined but is uniroot's upper limit. — Rui Barradas, Sep 13 '22 at 17:57
it is defined, as `upp_s <- 1/2 - 1e-20` (although this is a bad idea, `(1/2 - 1e-20) == 1/2` due to floating point error ...) — Ben Bolker, Sep 13 '22 at 18:00
`f3()` starts returning `NaN` somewhere between 0.497 and 0.498 — Ben Bolker, Sep 13 '22 at 18:04

Ben Bolker · Accepted Answer · 2022-09-13T21:48:54.457

I was able to get most of the way there with

upp_s <- 0.497
plot(Vectorize(f_est), from = 0.2, to = 100)

Not only is 1/2 - epsilon exactly equal to 1/2 for values of epsilon that are too small (due to floating point error), I found that f3() gives NaN for values >= 0.498. Setting upp_s to 0.497 worked OK.
plot() applied to a function calls curve(), which needs a function that can take a vector of x values.
The curve broke with "f() values at end points not of opposite sign" if I started the curve from 0.1; I didn't dig in further and try to diagnose what was going wrong.

PS. It is generally more numerically stable and efficient to do computations directly on the log scale where possible. In this case, that means using

f2 <- function(s) { (-3/2)*log(1-2*s) + (8*s)/(1-2*s)  }

instead of

f1 <- function(s) {
  (1 - 2*s)^(-3/2)*exp((8*s)/(1-2*s))
}
f2_orig <- function(s) {
  log(f1(s))
}
## check
all.equal(f2(0.25), f2_orig(0.25))  ## TRUE

Doing this and setting the lower bound of uniroot() to -500 lets us get pretty close to the zero boundary (although it looks both analytically and computationally as though the function diverges to -∞ as x goes to 0).

f3 <- Deriv(f2, 's')
upp_s <- 1/2 - 1e-10
lwr_a <- -500

f_est <- function(x) {
  f3a <- function(s) { f3(s = s) - x}
  s_ <- uniroot(f3a,
                lower = lwr_a,
                upper = upp_s)$root
  
  return(s_)
}
plot(Vectorize(f_est), from = 0.005, to = 100, log = "x")

You can also solve this analytically, or ask caracas (an R interface to sympy) to do it for you:

library(caracas)
x <- symbol("x"); s <- symbol("s")
## peek at f3() guts to find the expression for the derivative;
##  could also do the whole thing in caracas/sympy
solve_sys((11 +16*(s/(1-s*2)))/(1-s*2), x, list(s))
sol <- function(x) { (2*x - sqrt(32*x + 9) -3)/(4*x) }
curve(sol, add = TRUE, col = 2)

Why can't I use uniroot with plot?

1 Answers1