[Rd] uniroot violates bounds?
Martin Maechler
m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Wed Feb 22 11:32:25 CET 2023
>>>>> Martin Maechler on Mon, 20 Feb 2023 19:25:02 +0100 writes:
>>>>> Serguei Sokol via R-devel
>>>>> on Mon, 20 Feb 2023 10:59:22 +0100 writes:
>> Le 18/02/2023 à 21:44, J C Nash a écrit :
>>> I wrote first cut at unirootR for Martin M and he revised and put in
>>> Rmpfr.
>>>
>>> The following extends Ben's example, but adds the unirootR with trace
>>> output.
>>>
>>> c1 <- 4469.822
>>> c2 <- 572.3413
>>> f <- function(x) { c1/x - c2/(1-x) }; uniroot(f, c(1e-6, 1))
>>> uniroot(f, c(1e-6, 1))
>>> library(Rmpfr)
>>> unirootR(f, c(1e-6, 1), extendInt="no", trace=1)
>>>
>>> This gives more detail on the iterations, and it looks like the Inf is
>>> the
>>> issue. But certainly we could do more to avoid "gotchas" like this. If
>>> someone is interested in some back and forth, I'd be happy to give it a
>>> try, but I think progress would be better with more than one contributor.
>> For me, the following fix makes the job :
>> --- nlm.R.old 2018-09-25 10:44:49.000000000 +0200
>> +++ nlm.R 2023-02-20 10:46:39.893542531 +0100
>> @@ -143,14 +143,14 @@
>> if(check.conv) {
>> val <- tryCatch(.External2(C_zeroin2, function(arg) f(arg, ...),
>> - lower, upper, f.lower, f.upper,
>> + lower, upper, f.low., f.upp.,
>> tol, as.integer(maxiter)),
>> warning = function(w)w)
>> if(inherits(val, "warning"))
>> stop("convergence problem in zero finding: ",
>> conditionMessage(val))
>> } else {
>> val <- .External2(C_zeroin2, function(arg) f(arg, ...),
>> - lower, upper, f.lower, f.upper,
>> + lower, upper, f.low., f.upp.,
>> tol, as.integer(maxiter))
>> }
>> iter <- as.integer(val[2L])
>> Best,
>> Serguei.
> Thank you, Ben, John and Serguei !
> Serguei's fix shows the way -- even though it is not sufficient in general:
> f.low. and f.upp. may not be current after interval extensions
> in the doX case.
> But that's easy to amend, too.
> I will need to do a few more checks / code reading before I'm
> 100% convinced that such a change is good to be committed, but
> I'm pretty confident currently that it *is* ..
I have committed the changes to R-devel now
------------------------------------------------------------------------
r83889 | maechler | 2023-02-22 10:14:19 +0100
Changed paths:
M doc/NEWS.Rd
M src/library/stats/R/nlm.R
M src/library/stats/man/uniroot.Rd
M tests/reg-tests-1e.R
uniroot() close to |f()| = Inf staying in interval
------------------------------------------------------------------------
... in the NEWS thanking Ben and Serguei.
There's a small possibility that this did not only fix the bug,
but has also very slightly changed uniroot() results in
situations they were good already (so test code which over zealously
compares previous with current results may trigger ..).
but I think the extra safe guarding against +/- Inf values in
f(.) is well worth it
{The commit also improves the help page saying more about
'init.it' in the result.}
Martin
More information about the R-devel
mailing list