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 <profjcn...@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@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@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@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.