Uniroot solution in R

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?


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)

enter image description here

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)

enter image description here

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.

Try using a small interval but allow uniroot() to extend the interval:

uniroot(f, lower=0, upper=1, extendInt = "yes")$root
[1] -102.9519