[Rd] RFC: a "safe" uniroot() function for future R

Martin Maechler maechler at stat.math.ethz.ch
Thu May 30 12:49:33 CEST 2013

>>>>> Duncan Murdoch <murdoch.duncan at gmail.com>
>>>>>     on Thu, 30 May 2013 05:27:57 -0400 writes:

    > On Thu, May 30, 2013 at 4:18 AM, Martin Maechler <maechler at 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.

I did originally, too, but then became less sure about possible CRAN
checking "fallout"...

Merging the two functions into one, and not keeping the original
might be even the best solution, also easiest to maintain,
including documentation.

    >> 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:

Thanks a lot, Duncan!

    > 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?

Good point.  Note that it's completely unused if f() changes sign.
In that case and if  f() is "down crossing" at x0 , 
it would currently *not* warn even if  Sig == 1.
For me, using  Upcrossing = TRUE (or similar) would rather
suggest differently, i.e., I'd expect a warning when that
requirement is not fulfilled.

As you propose a default of NA, we *could* consider to warn in
such cases, where the user prescribes the slope-sign at the crossing...

    > 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?

good point.

    > 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.

yes, I agree.
And that's actually the point for discussion if we
should allow ourself to make this into uniroot():
The back compatibility is only guaranteed in those case where
old-uniroot() did *not* fail.

    > 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.

duh... of course!  
It shows I wrote the function for my own use, where example as the above
would not happen.

I've committed a simple change {and test} to catch such
To implement your proposal of *counting* the search
iterations among the total is trivial if we decide to change
the two functions into one uniroot(),
whereas with the current  unirootS() -> uniroot(), this would
need a bit of ugly code bloat... so I'd rather wait how this
RFC develops.


    > Duncan Murdoch


    >> -----------
    >> 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

More information about the R-devel mailing list