[Rd] uniroot violates bounds?

J C Nash pro|jcn@@h @end|ng |rom gm@||@com
Sat Feb 18 21:44:23 CET 2023


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.

Best,

John Nash

On 2023-02-18 12:28, Ben Bolker wrote:
> 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))
> 
> 
>     provides a root at -6.00e-05, which is outside of the specified bounds.  The default value of the "extendInt" 
> argument to uniroot() is "no", as far as I can see ...
> 
> $root
> [1] -6.003516e-05
> 
> $f.root
> [1] -74453981
> 
> $iter
> [1] 1
> 
> $init.it
> [1] NA
> 
> $estim.prec
> [1] 6.103516e-05
> 
> 
>    I suspect this fails because f(1) (value at the upper bound) is infinite, although setting interval to c(0.01, 1) 
> does work/give a sensible answer ...  (works for a lower bound of 1e-4, fails for 1e-5 ...)
> 
>    Setting the upper bound < 1 appears to avoid the problem.
> 
>   For what it's worth, the result has an "init.it" component, but the only thing the documentation says about it is " 
> component ‘init.it’ was added in R 3.1.0".
> 
>    And, I think (?) that the 'trace' argument only produces any output if the 'extendInt' option is enabled?
> 
>    Inspired by 
> https://stackoverflow.com/questions/75494696/solving-a-system-of-non-linear-equations-with-only-one-unknown/75494955#75494955
> 
>    cheers
>     Ben Bolker
> 
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



More information about the R-devel mailing list