With main R releases only happening yearly in spring, now is
good time to consider *and* discuss new features for what we
often call "R-devel" and more officially is
R Under development (unstable) (.....) -- "Unsuffered Consequences"
Here is one such example I hereby expose to public scrutiny:
A few minutes ago, I've committed the following to R-devel
(the 'trunk' in the svn repository for R):
------------------------------------------------------------------------
r62834 | maechler | 2013-05-30 10:01:33 +0200 (Thu, 30 May 2013) | 1 line
Changed paths:
M doc/NEWS.Rd
M src/library/stats/NAMESPACE
M src/library/stats/R/nlm.R
M src/library/stats/man/uniroot.Rd
M tests/Examples/stats-Ex.Rout.save
new "safe" uniroot() =: unirootS() [may change; see R-devel e-mail]
------------------------------------------------------------------------
The help file says
?unirootS()? is a ?safe? version of ?uniroot()?, built on
?uniroot()?, also useful as drop-in replacement of ?uniroot()? in
some cases. ?Safe? 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.
We've had this function, called safeUroot() in our package
copula for a while now, where an earlier not-exported version
has been in my package nor1mix even longer.
When I was tempted to introduce it into yet another CRAN package
I maintain, I decided that this may be a sign that such a
simple [ utility for / generalization of ] uniroot() should
probably rather be added to R itself.
The function definition, also visible in R's devel.sources, at the bottom of
https://svn.r-project.org/R/trunk/src/library/stats/R/nlm.R ,
shows that unirootS() is a wrapper for uniroot() and
is in principle 100% back compatible to uniroot() itself for all
the cases where f(lower) and f(upper) are of differing sign and
hence uniroot() does not give a quick error.
unirootS() just has three new optional arguments, all with their
defaults set such as to remain uniroot() compatible.
So, one option instead of the currently commited one would be to
adopt unirootS() as "the new uniroot()" and rename current
uniroot to .uniroot() {and still export both}.
The capital "S" in the function name and the 'Sig' name is of
course quite a matter of taste, and this case corresponds to my
taste, but even that is part of the RFC.
unirootS <- function(f, interval, ...,
lower = min(interval), upper = max(interval),
f.lower = f(lower, ...), f.upper = f(upper, ...),
Sig = NULL, check.conv = FALSE,
tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)
{
if ( is.null(Sig) && f.lower * f.upper > 0 ||
is.numeric(Sig) && (Sig*f.lower > 0 || Sig*f.upper < 0)) {
if(trace)
cat(sprintf("search in [%g,%g]%s", lower, upper,
if(trace >= 2)"\n" else " ... "))
Delta <- function(u) 0.01* pmax(1e-7, abs(u))
## Two cases:
if(is.null(Sig)) {
## case 1) 'Sig' unspecified --> extend (lower, upper) at the
same time
delta <- Delta(c(lower,upper))
while(isTRUE(f.lower*f.upper > 0) && any(iF <-
is.finite(c(lower,upper)))) {
if(iF[1]) f.lower <- f(lower <- lower - delta[1])
if(iF[2]) f.upper <- f(upper <- upper + delta[2])
if(trace >= 2)
cat(sprintf(" .. modified lower,upper: (%15g,%15g)\n",
lower,upper))
delta <- 2 * delta
}
} else {
## case 2) 'Sig' specified --> typically change only *one* of
lower, upper
## make sure we have Sig*f(lower) < 0 and Sig*f(upper) > 0:
delta <- Delta(lower)
while(isTRUE(Sig*f.lower > 0)) {
f.lower <- f(lower <- lower - delta)
if(trace) cat(sprintf(" .. modified lower: %g\n", lower))
delta <- 2 * delta
}
delta <- Delta(upper)
while(isTRUE(Sig*f.upper < 0)) {
f.upper <- f(upper <- upper + delta)
if(trace) cat(sprintf(" .. modified upper: %g\n", upper))
delta <- 2 * delta
}
}
if(trace && trace < 2)
cat(sprintf("extended to [%g, %g]\n", lower, upper))
}
if(!isTRUE(f.lower * f.upper <= 0))
stop("did not succeed extending the interval endpoints for f(lower) *
f(upper) <= 0")
if(check.conv) {
r <- tryCatch(uniroot(f, ..., lower=lower, upper=upper,
f.lower=f.lower, f.upper=f.upper,
tol=tol, maxiter=maxiter),
warning = function(w)w)
if(inherits(r, "warning"))
stop("convergence problem in zero finding: ", r$message)
else r
}
else
uniroot(f, ..., lower=lower, upper=upper,
f.lower=f.lower, f.upper=f.upper,
tol=tol, maxiter=maxiter)
}
-----------
As said, your comments are very welcome!
Note that I'm less interested in variations which gain 10-20% in
speed benchmarks, rather I'd appreciate proposals for changes
that give a "better" (in your sense) user interface.
Martin Maechler, ETH Zurich
On Thu, May 30, 2013 at 4:18 AM, Martin Maechler <maechler@stat.math.ethz.ch> wrote:> With main R releases only happening yearly in spring, now is > good time to consider *and* discuss new features for what we > often call "R-devel" and more officially is > R Under development (unstable) (.....) -- "Unsuffered Consequences" > > Here is one such example I hereby expose to public scrutiny: > > A few minutes ago, I've committed the following to R-devel > (the 'trunk' in the svn repository for R): > > ------------------------------------------------------------------------ > r62834 | maechler | 2013-05-30 10:01:33 +0200 (Thu, 30 May 2013) | 1 line > Changed paths: > M doc/NEWS.Rd > M src/library/stats/NAMESPACE > M src/library/stats/R/nlm.R > M src/library/stats/man/uniroot.Rd > M tests/Examples/stats-Ex.Rout.save > > new "safe" uniroot() =: unirootS() [may change; see R-devel e-mail] > ------------------------------------------------------------------------ > > The help file says > > ‘unirootS()’ is a “safe” version of ‘uniroot()’, built on > ‘uniroot()’, also useful as drop-in replacement of ‘uniroot()’ in > some cases. “Safe” 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. > > We've had this function, called safeUroot() in our package > copula for a while now, where an earlier not-exported version > has been in my package nor1mix even longer. > When I was tempted to introduce it into yet another CRAN package > I maintain, I decided that this may be a sign that such a > simple [ utility for / generalization of ] uniroot() should > probably rather be added to R itself. > > The function definition, also visible in R's devel.sources, at the bottom > of > https://svn.r-project.org/R/trunk/src/library/stats/R/nlm.R , > shows that unirootS() is a wrapper for uniroot() and > is in principle 100% back compatible to uniroot() itself for all > the cases where f(lower) and f(upper) are of differing sign and > hence uniroot() does not give a quick error. > unirootS() just has three new optional arguments, all with their > defaults set such as to remain uniroot() compatible. > > So, one option instead of the currently commited one would be to > adopt unirootS() as "the new uniroot()" and rename current > uniroot to .uniroot() {and still export both}. >I would probably prefer this.> > The capital "S" in the function name and the 'Sig' name is of > course quite a matter of taste, and this case corresponds to my > taste, but even that is part of the RFC. > > > unirootS <- function(f, interval, ..., > lower = min(interval), upper = max(interval), > f.lower = f(lower, ...), f.upper = f(upper, ...), > Sig = NULL, check.conv = FALSE, > tol = .Machine$double.eps^0.25, maxiter = 1000, trace > = 0) >A few comments: 1. I don't think the name "Sig" conveys the meaning of that parameter well. If specified, it determines whether the function is increasing or decreasing at the root, so maybe "Increasing" or "Upcrossing" (with a logical value, default NA) would be better? 2. In case 2 where the interval is expanded, wouldn't we save a bit of time by saving the initial values? E.g. if Sig == 1 so we want an upcrossing, but f.lower is positive, shouldn't we set upper to lower as we expand lower downwards? 3. Sometimes a user will want to force the solution to be between lower and upper, and will want to signal an error if they are not acceptable. If you do decide to merge this into uniroot that should be an option. 4. It should count the search for the interval among the iterations, and quit if it can't find an interval in that time. For example, unirootS( function(x) 1, c(0,1) ) never terminates. Duncan Murdoch {> if ( is.null(Sig) && f.lower * f.upper > 0 || > is.numeric(Sig) && (Sig*f.lower > 0 || Sig*f.upper < 0)) { > if(trace) > cat(sprintf("search in [%g,%g]%s", lower, upper, > if(trace >= 2)"\n" else " ... ")) > Delta <- function(u) 0.01* pmax(1e-7, abs(u)) > ## Two cases: > if(is.null(Sig)) { > ## case 1) 'Sig' unspecified --> extend (lower, upper) at the > same time > delta <- Delta(c(lower,upper)) > while(isTRUE(f.lower*f.upper > 0) && any(iF <- > is.finite(c(lower,upper)))) { > if(iF[1]) f.lower <- f(lower <- lower - delta[1]) > if(iF[2]) f.upper <- f(upper <- upper + delta[2]) > if(trace >= 2) > cat(sprintf(" .. modified lower,upper: (%15g,%15g)\n", > lower,upper)) > delta <- 2 * delta > } > } else { > ## case 2) 'Sig' specified --> typically change only *one* of > lower, upper > ## make sure we have Sig*f(lower) < 0 and Sig*f(upper) > 0: > delta <- Delta(lower) > while(isTRUE(Sig*f.lower > 0)) { > f.lower <- f(lower <- lower - delta) > if(trace) cat(sprintf(" .. modified lower: %g\n", lower)) > delta <- 2 * delta > } > delta <- Delta(upper) > while(isTRUE(Sig*f.upper < 0)) { > f.upper <- f(upper <- upper + delta) > if(trace) cat(sprintf(" .. modified upper: %g\n", upper)) > delta <- 2 * delta > } > } > if(trace && trace < 2) > cat(sprintf("extended to [%g, %g]\n", lower, upper)) > } > if(!isTRUE(f.lower * f.upper <= 0)) > stop("did not succeed extending the interval endpoints for > f(lower) * f(upper) <= 0") > if(check.conv) { > r <- tryCatch(uniroot(f, ..., lower=lower, upper=upper, > f.lower=f.lower, f.upper=f.upper, > tol=tol, maxiter=maxiter), > warning = function(w)w) > if(inherits(r, "warning")) > stop("convergence problem in zero finding: ", r$message) > else r > } > else > uniroot(f, ..., lower=lower, upper=upper, > f.lower=f.lower, f.upper=f.upper, > tol=tol, maxiter=maxiter) > } > > ----------- > > As said, your comments are very welcome! > Note that I'm less interested in variations which gain 10-20% in > speed benchmarks, rather I'd appreciate proposals for changes > that give a "better" (in your sense) user interface. > > Martin Maechler, ETH Zurich > > ______________________________________________ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >[[alternative HTML version deleted]]
Martin,
1. I'd vote for replacement.
2. The "Sig" argument was completely opaque to me. I'd vote
for something like this
bracket = "fixed" the old behavior, treat the upper and lower as
fixed values
= "search" expand the limits if needed
= "usearch" only the upper limit is expanded
= 'lsearch' only the lower
3. I agree with Duncan that the searches should count toward the total
iterations.
Terry T.
On 05/30/2013 05:00 AM, r-devel-request@r-project.org
wrote:> With main R releases only happening yearly in spring, now is
> good time to consider*and* discuss new features for what we
> often call "R-devel" and more officially is
> R Under development (unstable) (.....) -- "Unsuffered
Consequences"
>
> Here is one such example I hereby expose to public scrutiny:
>
> A few minutes ago, I've committed the following to R-devel
> (the 'trunk' in the svn repository for R):
>
> ------------------------------------------------------------------------
> r62834 | maechler | 2013-05-30 10:01:33 +0200 (Thu, 30 May 2013) | 1 line
> Changed paths:
> M doc/NEWS.Rd
> M src/library/stats/NAMESPACE
> M src/library/stats/R/nlm.R
> M src/library/stats/man/uniroot.Rd
> M tests/Examples/stats-Ex.Rout.save
>
> new "safe" uniroot() =: unirootS() [may change; see R-devel
e-mail]
> ------------------------------------------------------------------------
> ...
[[alternative HTML version deleted]]