10

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?

zx8754
  • 52,746
  • 12
  • 114
  • 209
shany
  • 175
  • 1
  • 1
  • 7
  • 5
    The error message seems self explanatory. f(0) and f(1000) have the same sign so uniroot refuses to search that interval for a zero. – Frank Aug 15 '16 at 18:58
  • 1
    You can plot this function over the range 0 to 1000. `curve(f(x), 0, 1000)`. Doesn't look like it crosses 0. – MrFlick Aug 15 '16 at 19:35

2 Answers2

28

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.
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • 2
    It's not true that uniroot uses bisection. It uses Brent's algorithm. In the R help file for uniroot, it states: "Based on ‘zeroin.c’ in http://www.netlib.org/c/brent.shar." Reading the documentation at netlib.org/c, we find the following information: "file brent.shar for Brent's univariate minimizer and zero finder. by Oleg Keselyov May 23, 1991 ref G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical computations." – inhuretnakht Mar 01 '20 at 15:48
5

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

uniroot(f, lower=0, upper=1, extendInt = "yes")$root
[1] -102.9519
Joe
  • 8,073
  • 1
  • 52
  • 58
John I.
  • 51
  • 1
  • 1
  • 3
    Just want to tell people that `extendInt` does not always work so don't believe it does the magic. Taking the quadratic polynomial example in my answer for example, adding `extendInt = "yes"` gives another "iteration reaches limits" error. Because no matter how much you extend the interval, a sign change never happens. – Zheyuan Li Aug 02 '18 at 01:37