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?
I was able to get most of the way there with
upp_s <- 0.497
plot(Vectorize(f_est), from = 0.2, to = 100)
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.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)