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)
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!
- Check the
f
above, is it really doing the right thing you want? I doubt something is wrong withg
; you might put brackets in the wrong place? - Once you get
f
correct, usecurve
to inspect a proper interval where there a root exists. Then useuniroot
.
Try using a small interval but allow uniroot() to extend the interval:
uniroot(f, lower=0, upper=1, extendInt = "yes")$root
[1] -102.9519