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 at 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 at 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 at 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.
These results reflect my experience with this sort of problem. A couple of comments: 1) optimx package has a multistart wrapper. I probably should have written one for nlsr. Maybe Bernard and I should work on that. The issues are largely to make things resistant to silly inputs, which even the smart users (you know, the ones looking back from the mirror) introduce. 2) Sometimes using the bounds constraint capability in nlsr can be helpful, e.g., to ensure the exponent parameters are kept apart, can be useful. 3) Combining with Duncan's suggestion of solving for the linear parameters also helps. All of the above can be sensitive to particular data. Best, JN On 2020-05-07 5:41 p.m., Bernard McGarvey wrote:> 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 at 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 at 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 at 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.
Dear all. Thank you for your answers. I will try Duncan's approach (if I could manage it). The issue is that first part of my data (actually temperature) up to certain time approximately follow one exponential. After that, another process prevails and the temperature increase starts to be "explosive". That is why I used these two exponentials. As I have many experiments I wanted to perform the fit programmatically. Which leads me to the approach that in each cycle I perform a plot which I visually inspect. If I consider the fit satisfactory I keep results. If not, I perform the fit with different starting values until it is OK. I am aware that it is not optimal but should be easiest. Thank you again. Best regards Petr> -----Original Message----- > From: J C Nash <profjcnash at gmail.com> > Sent: Friday, May 8, 2020 12:00 AM > To: Bernard McGarvey <mcgarvey.bernard at comcast.net>; PIKAL Petr > <petr.pikal at precheza.cz>; r-help <r-help at r-project.org> > Subject: Re: [R] unstable results of nlxb fit > > These results reflect my experience with this sort of problem. > > A couple of comments: > > 1) optimx package has a multistart wrapper. I probably should have written > one for nlsr. Maybe Bernard and I should work on that. The issues are > largely > to make things resistant to silly inputs, which even the smart users (you > know, > the ones looking back from the mirror) introduce. > > 2) Sometimes using the bounds constraint capability in nlsr can be helpful, > e.g., to ensure the exponent parameters are kept apart, can be useful. > > 3) Combining with Duncan's suggestion of solving for the linear parameters > also helps. > > All of the above can be sensitive to particular data. > > Best, JN > > On 2020-05-07 5:41 p.m., Bernard McGarvey wrote: > > 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 at 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 at 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 at 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.