[R] unstable results of nlxb fit
Bernard McGarvey
mcg@rvey@bern@rd @end|ng |rom comc@@t@net
Thu May 7 23:41:38 CEST 2020
John/Petr, I think there is an issue between a global optimum and local optima. I added a multistart loop around the code to see if I could find different solutions. Here is the code I added (I am not a great coder so please excuse any inefficiencies in this code segment):
# Multistart approach
NT <- 100
Results <- matrix(data=NA, nrow = NT, ncol=5, dimnames=list(NULL,c("SS", "A", "B", "a", "b")))
A1 <- runif(NT,0,100)
B1 <- runif(NT,0,100)
a1 <- runif(NT,0.0,0.1)
b1 <- runif(NT,0.0,0.1)
for (I in 1:NT) {
if (A1[I] > B1[I]) { # Ensure that the A'a are always the lower so that nlxb() always converge to the same values
A0 <- A1[I]
a0 <- a1[I]
A1[I] <- B1[I]
a1[I] <- b1[I]
B1[I] <- A0
b1[I] <- a0
}
fit <- nlxb(tsmes ~ A*exp(a*plast) + B* exp(b*plast), data=temp,
start=list(A=A1[I], B=B1[I], a=a1[I], b=b1[I]))
ccc <- coef(fit)
Results[I,1] <- fit$ssquares
Results[I,2] <- ccc[1]
Results[I,3] <- ccc[2]
Results[I,4] <- ccc[3]
Results[I,5] <- ccc[4]
}
Results
What I found is that the minimum SS generated at each trial had two distinct values, 417.8 and 3359.2. The A,B,a, and b values when the SS was 417.8 were all the same but I got different values for the case where the minimal SS was 3359.2. This indicates that the SS=417.8 may be the global minimum solution whereas the others are local optima. Here are the iteration results for a 100 trial multistart:
Results
SS A B a b
[1,] 3359.2 8.3546e+03 6.8321e+00 -1.988226 2.6139e-02
[2,] 3359.2 8.2865e+03 6.8321e+00 -5.201735 2.6139e-02
[3,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[4,] 3359.2 6.8321e+00 7.7888e+02 0.026139 -7.2812e-01
[5,] 3359.2 -3.9020e+01 4.5852e+01 0.026139 2.6139e-02
[6,] 3359.2 6.8321e+00 2.6310e+02 0.026139 -1.8116e+00
[7,] 3359.2 -2.1509e+01 2.8341e+01 0.026139 2.6139e-02
[8,] 3359.2 -3.8075e+01 4.4908e+01 0.026139 2.6139e-02
[9,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[10,] 3359.2 1.2466e+04 6.8321e+00 -4.196000 2.6139e-02
[11,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[12,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[13,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[14,] 3359.2 3.8018e+02 6.8321e+00 -0.806414 2.6139e-02
[15,] 3359.2 -3.1921e+00 1.0024e+01 0.026139 2.6139e-02
[16,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[17,] 3359.2 -1.5938e+01 2.2770e+01 0.026139 2.6139e-02
[18,] 3359.2 -3.1205e+01 3.8037e+01 0.026139 2.6139e-02
[19,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[20,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[21,] 3359.2 8.6627e+03 6.8321e+00 -3.319778 2.6139e-02
[22,] 3359.2 6.8321e+00 1.9318e+01 0.026139 -6.5773e-01
[23,] 3359.2 6.2991e+01 -5.6159e+01 0.026139 2.6139e-02
[24,] 3359.2 2.8865e-03 6.8321e+00 -1.576307 2.6139e-02
[25,] 3359.2 -1.2496e+01 1.9328e+01 0.026139 2.6139e-02
[26,] 3359.2 -5.9432e+00 1.2775e+01 0.026139 2.6139e-02
[27,] 3359.2 1.6884e+02 6.8321e+00 -211.866423 2.6139e-02
[28,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[29,] 3359.2 5.4972e+03 6.8321e+00 -3.432094 2.6139e-02
[30,] 3359.2 6.8321e+00 1.4427e+03 0.026139 -4.2771e+02
[31,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[32,] 3359.2 3.5760e+01 -2.8928e+01 0.026139 2.6139e-02
[33,] 3359.2 6.8321e+00 -4.0737e+02 0.026139 -6.7152e-01
[34,] 3359.2 6.8321e+00 1.2638e+04 0.026139 -2.8070e+00
[35,] 3359.2 1.1813e+01 -4.9807e+00 0.026139 2.6139e-02
[36,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[37,] 3359.2 6.8321e+00 1.2281e+03 0.026139 -3.0702e+02
[38,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[39,] 3359.2 -2.6850e+01 3.3682e+01 0.026139 2.6139e-02
[40,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[41,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[42,] 3359.2 -2.3279e+01 3.0111e+01 0.026139 2.6139e-02
[43,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[44,] 3359.2 6.8321e+00 1.4550e+03 0.026139 -4.0303e+00
[45,] 3359.2 -1.1386e+01 1.8218e+01 0.026139 2.6139e-02
[46,] 3359.2 8.8026e+02 6.8321e+00 -65.430608 2.6139e-02
[47,] 3359.2 -8.1985e+00 1.5031e+01 0.026139 2.6139e-02
[48,] 3359.2 -6.7767e+00 1.3609e+01 0.026139 2.6139e-02
[49,] 3359.2 -1.1436e+01 1.8268e+01 0.026139 2.6139e-02
[50,] 3359.2 1.0710e+04 6.8321e+00 -2.349659 2.6139e-02
[51,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[52,] 3359.2 6.8321e+00 7.1837e+02 0.026139 -7.4681e-01
[53,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[54,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[55,] 3359.2 -4.8774e+00 6.8321e+00 -16.405584 2.6139e-02
[56,] 3359.2 1.2687e+03 6.8321e+00 -3.775998 2.6139e-02
[57,] 3359.2 1.5529e+01 -8.6967e+00 0.026139 2.6139e-02
[58,] 3359.2 -1.0003e+01 1.6835e+01 0.026139 2.6139e-02
[59,] 3359.2 6.8321e+00 3.9291e+02 0.026139 -4.1974e+02
[60,] 3359.2 -2.1880e+01 2.8712e+01 0.026139 2.6139e-02
[61,] 3359.2 4.1736e+03 6.8321e+00 -10.711457 2.6139e-02
[62,] 3359.2 -3.3185e+01 4.0017e+01 0.026139 2.6139e-02
[63,] 3359.2 7.6732e+02 6.8321e+00 -0.723977 2.6139e-02
[64,] 3359.2 1.5334e+04 6.8321e+00 -52.573620 2.6139e-02
[65,] 3359.2 -2.9556e+01 3.6388e+01 0.026139 2.6139e-02
[66,] 3359.2 -1.0447e+00 7.8767e+00 0.026139 2.6139e-02
[67,] 3359.2 6.8321e+00 2.1471e+02 0.026139 -7.0582e+01
[68,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[69,] 3359.2 -2.2293e+01 2.9126e+01 0.026139 2.6139e-02
[70,] 3359.2 6.2259e+02 6.8321e+00 -2.782527 2.6139e-02
[71,] 3359.2 -1.4639e+01 2.1471e+01 0.026139 2.6139e-02
[72,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[73,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[74,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[75,] 3359.2 -2.3449e+01 3.0281e+01 0.026139 2.6139e-02
[76,] 3359.2 -2.5926e+01 6.8321e+00 -0.663656 2.6139e-02
[77,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[78,] 3359.2 6.8321e+00 6.9426e+02 0.026139 -1.9442e+00
[79,] 3359.2 2.8684e+02 6.8321e+00 -0.854394 2.6139e-02
[80,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[81,] 3359.2 -4.5066e+01 5.1899e+01 0.026139 2.6139e-02
[82,] 3359.2 4.4678e+03 6.8321e+00 -2.109446 2.6139e-02
[83,] 3359.2 3.1376e+03 6.8321e+00 -1.104803 2.6139e-02
[84,] 3359.2 6.8321e+00 1.1167e+02 0.026139 -1.0280e+00
[85,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[86,] 3359.2 5.3864e+02 6.8321e+00 -0.657971 2.6139e-02
[87,] 3359.2 4.8227e+01 6.8321e+00 -2.304024 2.6139e-02
[88,] 3359.2 -2.2048e+01 2.8880e+01 0.026139 2.6139e-02
[89,] 417.8 3.9452e-13 9.7727e+00 0.280227 2.1798e-02
[90,] 3359.2 6.8321e+00 -4.1689e+01 0.026139 -3.6049e+00
[91,] 417.8 9.7727e+00 3.9452e-13 0.021798 2.8023e-01
[92,] 3359.2 -4.1265e+01 4.8097e+01 0.026139 2.6139e-02
[93,] 3359.2 -1.1565e+01 1.8397e+01 0.026139 2.6139e-02
[94,] 3359.2 2.3698e+01 -1.6866e+01 0.026139 2.6139e-02
[95,] 3359.2 4.4700e+03 6.8321e+00 -12.836180 2.6139e-02
[96,] 3359.2 4.6052e+04 6.8321e+00 -7.158584 2.6139e-02
[97,] 3359.2 2.5464e+03 6.8321e+00 -1.811626 2.6139e-02
[98,] 3359.2 6.8321e+00 1.0338e+03 0.026139 -1.5365e+01
[99,] 3359.2 1.3783e+01 -6.9507e+00 0.026139 2.6139e-02
[100,] 3359.2 6.8321e+00 6.7153e+02 0.026139 -1.5975e+03
Hope this helps,
Bernard McGarvey
Director, Fort Myers Beach Lions Foundation, Inc.
Retired (Lilly Engineering Fellow).
> On May 7, 2020 at 9:33 AM J C Nash <profjcnash using gmail.com> wrote:
>
>
> 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.
> >
>
> ______________________________________________
> 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