[R] unstable results of nlxb fit

J C Nash pro|jcn@@h @end|ng |rom gm@||@com
Thu May 7 15:33:32 CEST 2020


The double exponential is well-known as a disaster to fit. Lanczos in his
1956 book Applied Analysis, p. 276 gives a good example which is worked through.
I've included it with scripts using nlxb in my 2014 book on Nonlinear Parameter
Optimization Using R Tools (Wiley). The scripts were on Wiley's site for the book,
but I've had difficulty getting Wiley to fix things and not checked lately if it
is still accessible. Ask off-list if you want the script and I'll dig into my
archives.

nlxb (preferably from nlsr which you used rather than nlmrt which is now not
maintained), will likely do as well as any general purpose code. There may be
special approaches that do a bit better, but I suspect the reality is that
the underlying problem is such that there are many sets of parameters with
widely different values that will get quite similar sums of squares.

Best, JN


On 2020-05-07 9:12 a.m., PIKAL Petr wrote:
> Dear all
> 
> I started to use nlxb instead of nls to get rid of singular gradient error.
> I try to fit double exponential function to my data, but results I obtain
> are strongly dependent on starting values. 
> 
> tsmes ~ A*exp(a*plast) + B* exp(b*plast)
> 
> Changing b from 0.1 to 0.01 gives me completely different results. I usually
> check result by a plot but could the result be inspected if it achieved good
> result without plotting?
> 
> Or is there any way how to perform such task?
> 
> Cheers
> Petr
> 
> Below is working example.
> 
>> dput(temp)
> temp <- structure(list(tsmes = c(31, 32, 32, 32, 32, 32, 32, 32, 33, 
> 34, 35, 35, 36, 36, 36, 37, 38, 39, 40, 40, 40, 40, 40, 41, 43, 
> 44, 44, 44, 46, 47, 47, 47, 47, 48, 49, 51, 51, 51, 52, 53, 54, 
> 54, 55, 57, 57, 57, 59, 59, 60, 62, 63, 64, 65, 66, 66, 67, 67, 
> 68, 70, 72, 74, 76, 78, 81, 84, 85, 86, 88, 90, 91, 92, 94, 96, 
> 97, 99, 100, 102, 104, 106, 109, 112, 115, 119, 123, 127, 133, 
> 141, 153, 163, 171), plast = c(50, 51, 52, 52, 53, 53, 53, 54, 
> 55, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 64, 64, 65, 65, 66, 
> 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 75, 76, 76, 77, 77, 78, 
> 78, 79, 80, 81, 82, 83, 84, 85, 85, 86, 86, 87, 88, 88, 89, 90, 
> 91, 91, 93, 93, 94, 95, 96, 96, 97, 98, 98, 99, 100, 100, 101, 
> 102, 103, 103, 104, 105, 106, 107, 107, 108, 109, 110, 111, 112, 
> 112, 113, 113, 114, 115, 116)), row.names = 2411:2500, class = "data.frame")
> 
> library(nlsr)
> 
> fit <- nlxb(tsmes ~ A*exp(a*plast) + B* exp(b*plast), data=temp,
> start=list(A=1, B=15, a=0.025, b=0.01))
> coef(fit)
>            A            B            a            b 
> 3.945167e-13 9.772749e+00 2.802274e-01 2.179781e-02 
> 
> plot(temp$plast, temp$tsmes, ylim=c(0,200))
> lines(temp$plast, predict(fit, newdata=temp), col="pink", lwd=3)
> ccc <- coef(fit)
> lines(0:120,ccc[1]*exp(ccc[3]*(0:120)))
> lines(0:120,ccc[2]*exp(ccc[4]*(0:120)), lty=3, lwd=2)
> 
> # wrong fit with slightly different b
> fit <- nlxb(tsmes ~ A*exp(a*plast) + B* exp(b*plast), data=temp,
> start=list(A=1, B=15, a=0.025, b=0.1))
> coef(fit)
>            A            B            a            b 
> 2911.6448377    6.8320597  -49.1373979    0.0261391 
> lines(temp$plast, predict(fit, newdata=temp), col="red", lwd=3)
> ccc <- coef(fit)
> lines(0:120,ccc[1]*exp(ccc[3]*(0:120)), col="red")
> lines(0:120,ccc[2]*exp(ccc[4]*(0:120)), lty=3, lwd=2, col="red")
> 
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



More information about the R-help mailing list