vuniroot: Vectorised One Dimensional Root (Zero) Finding
In rstpm2: Smooth Survival Models, Including Generalized Survival Models
vuniroot R Documentation
Vectorised One Dimensional Root (Zero) Finding
Description
The function vuniroot searches the interval from lower
to upper for a root (i.e., zero) of the vectorised function f with
respect to its first argument.
Setting extendInt to a non-"no" string, means searching
for the correct interval = c(lower,upper) if sign(f(x))
does not satisfy the requirements at the interval end points; see the
‘Details’ section.
Usage
vuniroot(f, interval, ...,
lower, upper,
f.lower = f(lower, ...), f.upper = f(upper, ...),
extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,
tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0,
n = NULL)
Arguments
f
the function for which the root is sought.
interval
a matrix with two columns containing the end-points of the interval
to be searched for the root.
...
additional named or unnamed arguments to be passed
to f
lower, upper
the lower and upper end points of the interval to
be searched.
f.lower, f.upper
the same as f(upper) and
f(lower), respectively. Passing these values from the caller
where they are often known is more economical as soon as f()
contains non-trivial computations.
extendInt
character string specifying if the interval
c(lower,upper) should be extended or directly produce an error
when f() does not have differing signs at the endpoints. The
default, "no", keeps the search interval and hence produces
an error. Can be abbreviated.
check.conv
logical indicating whether a convergence warning of the
underlying vuniroot should be caught as an error and if
non-convergence in maxiter iterations should be an error
instead of a warning.
tol
the desired accuracy (convergence tolerance).
maxiter
the maximum number of iterations.
trace
integer number; if positive, tracing information is
produced. Higher values giving more details.
n
integer number; size of input vector to f
(only used if lower and upper are of length 1)
Details
Note that arguments after ... must be matched exactly.
Either interval or both lower and upper must be
specified: the upper endpoint must be strictly larger than the lower
endpoint.
The function values at the endpoints must be of opposite signs (or
zero), for extendInt="no", the default. Otherwise, if
extendInt="yes", the interval is extended on both sides, in
search of a sign change, i.e., until the search interval [l,u]
satisfies f(l) \cdot f(u) \le 0.
If it is known how f changes sign at the root
x_0, that is, if the function is increasing or decreasing there,
extendInt can (and typically should) be specified as
"upX" (for “upward crossing”) or "downX",
respectively. Equivalently, define S := \pm 1, to
require S = \mathrm{sign}(f(x_0 + \epsilon)) at the solution. In that case, the search interval [l,u]
possibly is extended to be such that S\cdot f(l)\le 0 and S \cdot f(u) \ge 0.
vuniroot() uses a C++ subroutine based on ‘"zeroin"’ (from Netlib)
and algorithms given in the reference below. They assume a
continuous function (which then is known to have at least one root in
the interval).
Convergence is declared either if f(x) == 0 or the change in
x for one step of the algorithm is less than tol (plus an
allowance for representation error in x).
If the algorithm does not converge in maxiter steps, a warning
is printed and the current approximation is returned.
f will be called as f(x, ...) for a numeric value
of x.
The argument passed to f has special semantics and used to be
shared between calls. The function should not copy it.
Value
A list with at least three components: root and f.root
give the location of the root and the value of the function evaluated
at that point. iter gives the number of
iterations used.
Further components may be added in future: component init.it
was added in R 3.1.0.
Source
Based on ‘zeroin.c’ in https://netlib.org/c/brent.shar.
References
Brent, R. (1973)
Algorithms for Minimization without Derivatives.
Englewood Cliffs, NJ: Prentice-Hall.
See Also
uniroot for the standard single root solver
polyroot for all complex roots of a polynomial;
optimize, nlm.
Examples
require(utils) # for str
## some platforms hit zero exactly on the first step:
## if so the estimated precision is 2/3.
f <- function (x, a) x - a
str(xmin <- vuniroot(f, lower=c(0, 0), upper=c(1,1), tol = 0.0001, a = c(1/3,2/3)))
## same example with scalars for lower and upper -- using the n argument
str(xmin <- vuniroot(f, lower=0, upper=1, tol = 0.0001, n=2, a = c(1/3,2/3)))
## handheld calculator example: fixed point of cos(.):
vuniroot(function(x) cos(x) - x, lower = -pi, upper = pi, tol = 1e-9)$root
str(vuniroot(function(x) x*(x^2-1) + .5, lower = -2, upper = 2,
tol = 0.0001))
str(vuniroot(function(x) x*(x^2-1) + .5, lower = -2, upper = 2,
tol = 1e-10))
## Find the smallest value x for which exp(x) > 0 (numerically):
r <- vuniroot(function(x) 1e80*exp(x) - 1e-300, cbind(-1000, 0), tol = 1e-15)
str(r, digits.d = 15) # around -745, depending on the platform.
exp(r$root) # = 0, but not for r$root * 0.999...
minexp <- r$root * (1 - 10*.Machine$double.eps)
exp(minexp) # typically denormalized
##--- vuniroot() with new interval extension + checking features: --------------
f1 <- function(x) (121 - x^2)/(x^2+1)
f2 <- function(x) exp(-x)*(x - 12)
tools::assertCondition(vuniroot(f1, cbind(0,10)),
"error", verbose=TRUE)
tools::assertCondition(vuniroot(f2, cbind(0, 2)),
"error", verbose=TRUE)
##--> error: f() .. end points not of opposite sign
## where as 'extendInt="yes"' simply first enlarges the search interval:
u1 <- vuniroot(f1, cbind(0,10),extendInt="yes", trace=1)
u2 <- vuniroot(f2, cbind(0,2), extendInt="yes", trace=2)
stopifnot(all.equal(u1$root, 11, tolerance = 1e-5),
all.equal(u2$root, 12, tolerance = 6e-6))
## The *danger* of interval extension:
## No way to find a zero of a positive function, but
## numerically, f(-|M|) becomes zero :
tools::assertCondition(u3 <- vuniroot(exp, cbind(0,2), extendInt="yes", trace=TRUE),
"error", verbose=TRUE)
## Nonsense example (must give an error):
tools::assertCondition( vuniroot(function(x) 1, cbind(0,1), extendInt="yes"),
"error", verbose=TRUE)
## Convergence checking :
sinc_ <- function(x) ifelse(x == 0, 1, sin(x)/x)
curve(sinc_, -6,18); abline(h=0,v=0, lty=3, col=adjustcolor("gray", 0.8))
vuniroot(sinc_, cbind(0,5), extendInt="yes", maxiter=4) #-> "just" a warning
## now with check.conv=TRUE, must signal a convergence error :
vuniroot(sinc_, cbind(0,5), extendInt="yes", maxiter=4, check.conv=TRUE)
### Weibull cumulative hazard (example origin, Ravi Varadhan):
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u
n <- 10
u <- -log(runif(n))
a <- 1/2
b <- 1
## Find failure times
ru <- vuniroot(froot, u=u, a=a, b=b, interval= cbind(rep(1.e-14,n), rep(1e4,n)),
extendInt="yes")$root
ru2 <- vuniroot(froot, u=u, a=a, b=b, interval= cbind(rep(0.01,n), rep(10,n)),
extendInt="yes")$root
stopifnot(all.equal(ru, ru2, tolerance = 6e-6))
r1 <- vuniroot(froot, u= 0.99, a=a, b=b, interval= cbind(0.01, 10),
extendInt="up")
stopifnot(all.equal(0.99, cumhaz(r1$root, a=a, b=b)))
## An error if 'extendInt' assumes "wrong zero-crossing direction":
vuniroot(froot, u= 0.99, a=a, b=b, interval= cbind(0.1, 10), extendInt="down")
rstpm2 documentation built on Oct. 29, 2024, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| vuniroot | R Documentation |
Vectorised One Dimensional Root (Zero) Finding
Description
The function vuniroot searches the interval from lower
to upper for a root (i.e., zero) of the vectorised function f with
respect to its first argument.
Setting extendInt to a non-"no" string, means searching
for the correct interval = c(lower,upper) if sign(f(x))
does not satisfy the requirements at the interval end points; see the
‘Details’ section.
Usage
vuniroot(f, interval, ...,
lower, upper,
f.lower = f(lower, ...), f.upper = f(upper, ...),
extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,
tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0,
n = NULL)
Arguments
f |
the function for which the root is sought. |
interval |
a matrix with two columns containing the end-points of the interval to be searched for the root. |
... |
additional named or unnamed arguments to be passed
to |
lower, upper |
the lower and upper end points of the interval to be searched. |
f.lower, f.upper |
the same as |
extendInt |
character string specifying if the interval
|
check.conv |
logical indicating whether a convergence warning of the
underlying |
tol |
the desired accuracy (convergence tolerance). |
maxiter |
the maximum number of iterations. |
trace |
integer number; if positive, tracing information is produced. Higher values giving more details. |
n |
integer number; size of input vector to |
Details
Note that arguments after ... must be matched exactly.
Either interval or both lower and upper must be
specified: the upper endpoint must be strictly larger than the lower
endpoint.
The function values at the endpoints must be of opposite signs (or
zero), for extendInt="no", the default. Otherwise, if
extendInt="yes", the interval is extended on both sides, in
search of a sign change, i.e., until the search interval [l,u]
satisfies f(l) \cdot f(u) \le 0.
If it is known how f changes sign at the root
x_0, that is, if the function is increasing or decreasing there,
extendInt can (and typically should) be specified as
"upX" (for “upward crossing”) or "downX",
respectively. Equivalently, define S := \pm 1, to
require S = \mathrm{sign}(f(x_0 + \epsilon)) at the solution. In that case, the search interval [l,u]
possibly is extended to be such that S\cdot f(l)\le 0 and S \cdot f(u) \ge 0.
vuniroot() uses a C++ subroutine based on ‘"zeroin"’ (from Netlib)
and algorithms given in the reference below. They assume a
continuous function (which then is known to have at least one root in
the interval).
Convergence is declared either if f(x) == 0 or the change in
x for one step of the algorithm is less than tol (plus an
allowance for representation error in x).
If the algorithm does not converge in maxiter steps, a warning
is printed and the current approximation is returned.
f will be called as f(x, ...) for a numeric value
of x.
The argument passed to f has special semantics and used to be
shared between calls. The function should not copy it.
Value
A list with at least three components: root and f.root
give the location of the root and the value of the function evaluated
at that point. iter gives the number of
iterations used.
Further components may be added in future: component init.it
was added in R 3.1.0.
Source
Based on ‘zeroin.c’ in https://netlib.org/c/brent.shar.
References
Brent, R. (1973) Algorithms for Minimization without Derivatives. Englewood Cliffs, NJ: Prentice-Hall.
See Also
uniroot for the standard single root solver
polyroot for all complex roots of a polynomial;
optimize, nlm.
Examples
require(utils) # for str
## some platforms hit zero exactly on the first step:
## if so the estimated precision is 2/3.
f <- function (x, a) x - a
str(xmin <- vuniroot(f, lower=c(0, 0), upper=c(1,1), tol = 0.0001, a = c(1/3,2/3)))
## same example with scalars for lower and upper -- using the n argument
str(xmin <- vuniroot(f, lower=0, upper=1, tol = 0.0001, n=2, a = c(1/3,2/3)))
## handheld calculator example: fixed point of cos(.):
vuniroot(function(x) cos(x) - x, lower = -pi, upper = pi, tol = 1e-9)$root
str(vuniroot(function(x) x*(x^2-1) + .5, lower = -2, upper = 2,
tol = 0.0001))
str(vuniroot(function(x) x*(x^2-1) + .5, lower = -2, upper = 2,
tol = 1e-10))
## Find the smallest value x for which exp(x) > 0 (numerically):
r <- vuniroot(function(x) 1e80*exp(x) - 1e-300, cbind(-1000, 0), tol = 1e-15)
str(r, digits.d = 15) # around -745, depending on the platform.
exp(r$root) # = 0, but not for r$root * 0.999...
minexp <- r$root * (1 - 10*.Machine$double.eps)
exp(minexp) # typically denormalized
##--- vuniroot() with new interval extension + checking features: --------------
f1 <- function(x) (121 - x^2)/(x^2+1)
f2 <- function(x) exp(-x)*(x - 12)
tools::assertCondition(vuniroot(f1, cbind(0,10)),
"error", verbose=TRUE)
tools::assertCondition(vuniroot(f2, cbind(0, 2)),
"error", verbose=TRUE)
##--> error: f() .. end points not of opposite sign
## where as 'extendInt="yes"' simply first enlarges the search interval:
u1 <- vuniroot(f1, cbind(0,10),extendInt="yes", trace=1)
u2 <- vuniroot(f2, cbind(0,2), extendInt="yes", trace=2)
stopifnot(all.equal(u1$root, 11, tolerance = 1e-5),
all.equal(u2$root, 12, tolerance = 6e-6))
## The *danger* of interval extension:
## No way to find a zero of a positive function, but
## numerically, f(-|M|) becomes zero :
tools::assertCondition(u3 <- vuniroot(exp, cbind(0,2), extendInt="yes", trace=TRUE),
"error", verbose=TRUE)
## Nonsense example (must give an error):
tools::assertCondition( vuniroot(function(x) 1, cbind(0,1), extendInt="yes"),
"error", verbose=TRUE)
## Convergence checking :
sinc_ <- function(x) ifelse(x == 0, 1, sin(x)/x)
curve(sinc_, -6,18); abline(h=0,v=0, lty=3, col=adjustcolor("gray", 0.8))
vuniroot(sinc_, cbind(0,5), extendInt="yes", maxiter=4) #-> "just" a warning
## now with check.conv=TRUE, must signal a convergence error :
vuniroot(sinc_, cbind(0,5), extendInt="yes", maxiter=4, check.conv=TRUE)
### Weibull cumulative hazard (example origin, Ravi Varadhan):
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u
n <- 10
u <- -log(runif(n))
a <- 1/2
b <- 1
## Find failure times
ru <- vuniroot(froot, u=u, a=a, b=b, interval= cbind(rep(1.e-14,n), rep(1e4,n)),
extendInt="yes")$root
ru2 <- vuniroot(froot, u=u, a=a, b=b, interval= cbind(rep(0.01,n), rep(10,n)),
extendInt="yes")$root
stopifnot(all.equal(ru, ru2, tolerance = 6e-6))
r1 <- vuniroot(froot, u= 0.99, a=a, b=b, interval= cbind(0.01, 10),
extendInt="up")
stopifnot(all.equal(0.99, cumhaz(r1$root, a=a, b=b)))
## An error if 'extendInt' assumes "wrong zero-crossing direction":
vuniroot(froot, u= 0.99, a=a, b=b, interval= cbind(0.1, 10), extendInt="down")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.