>>>>> J C Nash
>>>>> on Fri, 20 Aug 2021 11:06:25 -0400 writes:
> In our work on a Google Summer of Code project
> "Improvements to nls()", the code has proved sufficiently
> entangled that we have found (so far!) few
> straightforward changes that would not break legacy
> behaviour. One issue that might be fixable is that nls()
> returns no result if it encounters some computational
> blockage AFTER it has already found a much better "fit"
> i.e. set of parameters with smaller sum of squares. Here
> is a version of the Tetra example:
time=c( 1, 2, 3, 4, 6 , 8, 10, 12, 16)
conc = c( 0.7, 1.2, 1.4, 1.4, 1.1, 0.8, 0.6, 0.5, 0.3)
NLSdata <- data.frame(time,conc)
NLSstart <-c(lrc1=-2,lrc2=0.25,A1=150,A2=50) # a starting vector (named!)
NLSformula <-conc ~ A1*exp(-exp(lrc1)*time)+A2*exp(-exp(lrc2)*time)
tryit <- try(nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE))
print(tryit)
> If you run this, tryit does not give information that the
> sum of squares has been reduced from > 60000 to < 2, as
> the trace shows.
> Should we propose that this be changed so the returned
> object gives the best fit so far, albeit with some form of
> message or return code to indicate that this is not
> necessarily a conventional solution? Our concern is that
> some examples might need to be adjusted slightly, or we
> might simply add the "try-error" class to the output
> information in such cases.
> Comments are welcome, as this is as much an infrastructure
> matter as a computational one.
Hmm... many years ago, we had introduced the 'warnOnly=TRUE'
option to nls() i.e., nls.control() exactly for such cases,
where people would still like to see the solution:
So,
------------------------------------------------------------------------------
> try2 <- nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE,
control = nls.control(warnOnly=TRUE))
61215.76 (3.56e+03): par = (-2 0.25 150 50)
2.175672 (2.23e+01): par = (-1.9991 0.3171134 2.618224 -1.366768)
1.621050 (7.14e+00): par = (-1.960475 -2.620293 2.575261 -0.5559918)
Warning message:
In nls(NLSformula, data = NLSdata, start = NLSstart, trace = TRUE, :
singular gradient
> try2
Nonlinear regression model
model: conc ~ A1 * exp(-exp(lrc1) * time) + A2 * exp(-exp(lrc2) * time)
data: NLSdata
lrc1 lrc2 A1 A2
-22.89 96.43 156.70 -156.68
residual sum-of-squares: 218483
Number of iterations till stop: 2
Achieved convergence tolerance: 7.138
Reason stopped: singular gradient
> coef(try2)
lrc1 lrc2 A1 A2
-22.88540 96.42686 156.69547 -156.68461
> summary(try2)
Error in chol2inv(object$m$Rmat()) :
element (3, 3) is zero, so the inverse cannot be computed
>
------------------------------------------------------------------------------
and similar for vcov(), of course, where the above error
originates.
{ I think GSoC (andr other) students should start by studying and
exploring relevant help pages before drawing conclusions
......
but yes, I've been born in the last millennium ...
}
;-)
Have a nice weekend!
Martin
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
