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.
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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
