# Issue

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?

# Solution

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)
```

Answered By - Ben Bolker Answer Checked By - David Marino (PHPFixing Volunteer)

## 0 Comments:

## Post a Comment

Note: Only a member of this blog may post a comment.