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
RFC: a "safe" uniroot() function for future R
7 messages · Duncan Murdoch, Martin Maechler, Ravi Varadhan
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20130530/d70adccf/attachment.pl>
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
examples.
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.
Martin
> 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
Dear Martin, I am not sure I like this idea of expanding the interval. It can have bad consequences. The best feature of uniroot is that it makes the user think about the behavior of the function. Your suggestion is in the spirit of making him unthink (if there is such a word!). Here is a cautionary example:
f <- function(x) exp(-x)
unirootS(f, c(0,2))
$root [1] 1312.7 $f.root [1] 0 $iter [1] 0 $estim.prec [1] 0 The existing `uniroot' does the right thing.
uniroot(f, c(0,2))
Error in uniroot(f, c(0, 2)) : f() values at end points not of opposite sign
Best, Ravi
-----Original Message----- From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org] On Behalf Of Duncan Murdoch Sent: Thursday, May 30, 2013 5:28 AM To: Martin Maechler Cc: R. Devel List Subject: Re: [Rd] RFC: a "safe" uniroot() function for future R 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.
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 at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[[alternative HTML version deleted]]
On 13-05-30 6:49 AM, Martin Maechler wrote:
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.
I don't think that's true (but maybe you've changed it since I downloaded; I'm offline right now). For example, try unirootS(function(x) x, c(-1, 1), Sig=-1) This fails after 1000 iterations, whereas Sig=+1 is fine. Duncan Murdoch
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
examples.
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.
Martin
> 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
Thank you, Ravi and Therry,
Ravi Varadhan <ravi.varadhan at jhu.edu>
on Thu, 30 May 2013 14:20:19 +0000 writes:
> Dear Martin,
> I am not sure I like this idea of expanding the interval. It can have bad consequences. The best feature of uniroot is that it makes the user think about the behavior of the function. Your suggestion is in the spirit of making him unthink (if there is such a word!).
Not necessarily:
My use case is inversion of a function that I know to be
monotone, and I know F^{-1}(a) has one well defined solution,
but sometimes my guess for an intervall is only approximate and
then simple uniroot() will fail. I've seen many such cases.
So, to find F^{-1}(a)
I really want a version of uniroot()
uniroot(function(x) F(x) - a, *)
where my initial interval is only approximate and may be too small.
> Here is a cautionary example:
>> f <- function(x) exp(-x)
>> unirootS(f, c(0,2))
> $root
> [1] 1312.7
> $f.root
> [1] 0
> $iter
> [1] 0
> $estim.prec
> [1] 0
yes; even shorter and equivalent is
unirootS(exp, c(-2,0))
> The existing `uniroot' does the right thing.
>> uniroot(f, c(0,2))
> Error in uniroot(f, c(0, 2)) :
> f() values at end points not of opposite sign
yes... but
uniroot(exp, c(-750, 0))
also gives a result like the above.
As Duncan said, and I agree, we will keep the option to say
"do use the provided interval and do not enlarge it".
With the current version, this would be with 'Sig = 0' :
unirootS(exp, c(-2, 0), Sig=0)
Error in unirootS(exp, c(-2, 0), Sig = 0) :
f() values at end points not of opposite sign
-------
We may think of extra arguments to limit the "search outside" in
some way, e.g., by an extra 'maxit' for these steps that could
be set low.
In that .. and maybe in any case, we probably should also
consider changing the way the "initial step size" delta is computed.
I now think that the initial delta should be something like
(upper - lower) / 16
rather than the current way, where the delta() is
computed separately and independently for 'lower' and 'upper'.
Martin
> Best,
> Ravi
>> -----Original Message-----
>> From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org]
>> On Behalf Of Duncan Murdoch
>> Sent: Thursday, May 30, 2013 5:28 AM
>> To: Martin Maechler
>> Cc: R. Devel List
>> Subject: Re: [Rd] RFC: a "safe" uniroot() function for future R
>>
>> 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.
>>
>>
>> >
>> > 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
On 13-05-30 6:49 AM, Martin Maechler wrote:
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.
I don't think that's true (but maybe you've changed it since I downloaded; I'm offline right now). For example, try
unirootS(function(x) x, c(-1, 1), Sig=-1)
This fails after 1000 iterations, whereas Sig=+1 is fine.
Indeed. You are entirely right, and my note was wrong.
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.
so the above is moot.
From the current feedbacks, I'd come to propose / further
discuss the following issues:
1) the goal is to remain with one function uniroot()
2) Instead of the 'Sig' = "sign(f'(x_0))" {not quite, but typically}
with 4 different value classes, namely
NULL, -1, 0, 1, (where +/- 1 are equivalent to any positive
or negative finite number respectively),
we should either use a string with 4 different possible values
or a {logical or NULL}, say 'upcrossing'
(which also gives 4 values, NULL, NA, FALSE, TRUE).
3) [I'm not sure about this:]
The new default of the 'Sig' replacement would correspond to
the current Sig = NULL which does extend the search
interval when that does not constitute a sign change.
Alternatively, implicitly proposed by Ravi Varadhan, the
default would correspond to Sig = 0, i.e. to the current
uniroot() behavior which signals an error as soon as the initial
interval is not large enough.
Further opinions and suggestions for '2)' and '3)' are still
very welcome!
Martin
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
examples.
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.
Martin
> 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