Uniroot solution in R

Question

I would like to find the root of the following function:

       x=0.5
       f <- function(y) ((1-pbeta(1-exp(-0.002926543
        *( 107.2592+y)^1.082618 *exp(0.04097536*(107.2592+y))),shape1=0.2640229,shape2=0.1595841)) -
(1-pbeta(1-exp(-0.002926543*(x)^1.082618 *exp(0.04097536*(x))),shape1=0.2640229,shape2=0.1595841))^2)

sroot=uniroot(f, lower=0, upper=1000)$root

Error in uniroot(f, lower = 0, upper = 1000) : f() values at end points not of opposite sign

How can I solve the error?

Answer 1

uniroot() and caution of its use

uniroot is implementing the crude bisection method. Such method is much simpler that (quasi) Newton's method, but need stronger assumption to ensure the existence of a root: f(lower) * f(upper) < 0.

This can be quite a pain,as such assumption is a sufficient condition, but not a necessary one. In practice, if f(lower) * f(upper) > 0, it is still possible that a root exist, but since this is not of 100% percent sure, bisection method can not take the risk.

Consider this example:

# a quadratic polynomial with root: -2 and 2
f <- function (x) x ^ 2 - 4

Obviously, there are roots on [-5, 5]. But

uniroot(f, lower = -5, upper = 5)
#Error in uniroot(f, lower = -5, upper = 5) : 
#  f() values at end points not of opposite sign

In reality, the use of bisection method requires observation / inspection of f, so that one can propose a reasonable interval where root lies. In R, we can use curve():

curve(f, from = -5, to = 5); abline(h = 0, lty = 3)

From the plot, we observe that a root exist in [-5, 0] or [0, 5]. So these work fine:

uniroot(f, lower = -5, upper = 0)
uniroot(f, lower = 0, upper = 5)

---

Your problem

Now let's try your function (I have split it into several lines for readability; it is also easy to check correctness this way):

f <- function(y) {
  g <- function (u) 1 - exp(-0.002926543 * u^1.082618 * exp(0.04097536 * u))
  a <- 1 - pbeta(g(107.2592+y), 0.2640229, 0.1595841)
  b <- 1 - pbeta(g(x), 0.2640229, 0.1595841)
  a - b^2
  }

x <- 0.5
curve(f, from = 0, to = 1000)

How could this function be a horizontal line? It can't have a root!

1. Check the f above, is it really doing the right thing you want? I doubt something is wrong with g; you might put brackets in the wrong place?
2. Once you get f correct, use curve to inspect a proper interval where there a root exists. Then use uniroot.