# Starting estimates for nls Exponential Fit

11 messages
Open this post in threaded view
|

## Starting estimates for nls Exponential Fit

 Hello everyone, I have come across a bit of an odd problem: I am currently analysing data PCR reaction data of which part is behaving exponential. I would like to fit the exponential part to the following: y0 + alpha * E^t In which Y0 is the groundphase, alpha a fluorescence factor, E the efficiency of the reaction & t is time (in cycles) I can get this to work for most of my reactions, but part fails the nls fit (Convergence failure: iteration limit reached without convergence). This mainly has to do with the starting estimates for the nls function. I have used various 'smart' ways of calculating starting estimates (grid searches, optim(), etc.) but however close the starting estimates are to the actual values, nls still fails for many datasets. The weirdest thing is, I have one set of starting estimates (alpha= 0.5 and E=1.85) which are totaly arbitray and these DO work for, say 99% of the datasets. I don't know why this is and I don't why my 'good estimates' do not work. I am desperatly seeking a way of calculating working starting estimates for my nls function. Can anybody give me a hand? thanks, Anto R version 2.9.2 example dataset: ExponValues [1] 2018.34 2012.54 2018.85 2023.52 2054.58 2132.61 2247.17 2468.32 2778.47 ExponCycles [1] 17 18 19 20 21 22 23 24 25 Example starting estimate calculation Y0 <- ExponValues[1] E <- weighted.mean((ExponValues-eY0)[-1][-1]/(ExponValues-eY0)[-1][-(length(ExponValues)-1)],(1-abs(seq(-1,1,length=(length(ExponValues)-2)))^3)^3) alpha <-  weighted.mean((ExponValues[-1]-ExponValues[-length(ExponValues)])/((E^ExponCycles[-1])-(E^ExponCycles[-length(ExponCycles)])),(1-abs(seq(-1,1,length=(length(ExponCycles)-1)))^3)^3) Optional parameter optimisation: Esp <- optim(c(Y0=eY0,a=alpha,E=E),function(x){return(sum(abs(ExponValues-(x[1]+x[2]*x[3]^ExponCycles))))})\$par nls function: Emodel<-try(nls(ExponValues ~ (Y0 + a*(E^ExponCycles)) , start= Esp, algorithm="port"),silent=TRUE)
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 If you could reformulate your model as alpha * (y0 + E^t) then you could use nls with alg="plinear" (alpha then would be eliminated from the nonlinear param and treated as conditionally linear) and this would help with convergence. Else you can try package DEoptim to get the starting values; the advantage over optim is that you need then lower and upper bounds on the parameters, not more starting values; the bounds however should be appropriate and as tight as possible. Also: by default nls uses max. 50 iter.  Depending on where you start off you may need more than this. ## ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,                  2468.32,2778.47) ExponCycles <- c(17,18,19,20,21,22,23,24,25) library(DEoptim) mod <- function(x) x[1] + x[2]*x[3]^ExponCycles fun <- function(x) sum((ExponValues-mod(x))^2) ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),               control=list(trace=FALSE)) pa <- ss\$optim\$bestmem ## now have par est that give OK fit matplot(cbind(ExponValues, mod(pa)),type="l") names(pa) <- c("Y0", "a", "E") ## fit w/port and GN Emodel<- nls(ExponValues ~ (Y0 + a*(E^ExponCycles)),              start=pa, algorithm="port",              control=list(maxiter=1000, warnOnly=TRUE)) Emodel1 <- nls(ExponValues ~ (Y0 + a*(E^ExponCycles)), start=pa,              control=list(maxiter=1000, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodel), fitted(Emodel1)),type="l") On Tue, 1 Dec 2009, Anto wrote: > > Hello everyone, > > I have come across a bit of an odd problem: I am currently analysing data > PCR reaction data of which part is behaving exponential. I would like to fit > the exponential part to the following: > > y0 + alpha * E^t > > In which Y0 is the groundphase, alpha a fluorescence factor, E the > efficiency of the reaction & t is time (in cycles) > > I can get this to work for most of my reactions, but part fails the nls fit > (Convergence failure: iteration limit reached without convergence). This > mainly has to do with the starting estimates for the nls function. I have > used various 'smart' ways of calculating starting estimates (grid searches, > optim(), etc.) but however close the starting estimates are to the actual > values, nls still fails for many datasets. > > The weirdest thing is, I have one set of starting estimates (alpha= 0.5 and > E=1.85) which are totaly arbitray and these DO work for, say 99% of the > datasets. I don't know why this is and I don't why my 'good estimates' do > not work. I am desperatly seeking a way of calculating working starting > estimates for my nls function. Can anybody give me a hand? > > thanks, > > Anto > > R version 2.9.2 > > example dataset: > > ExponValues > [1] 2018.34 2012.54 2018.85 2023.52 2054.58 2132.61 2247.17 2468.32 2778.47 > > ExponCycles > [1] 17 18 19 20 21 22 23 24 25 > > > Example starting estimate calculation > > Y0 <- ExponValues[1] > E <- > weighted.mean((ExponValues-eY0)[-1][-1]/(ExponValues-eY0)[-1][-(length(ExponValues)-1)],(1-abs(seq(-1,1,length=(length(ExponValues)-2)))^3)^3) > alpha <- > weighted.mean((ExponValues[-1]-ExponValues[-length(ExponValues)])/((E^ExponCycles[-1])-(E^ExponCycles[-length(ExponCycles)])),(1-abs(seq(-1,1,length=(length(ExponCycles)-1)))^3)^3) > > > Optional parameter optimisation: > > Esp <- > optim(c(Y0=eY0,a=alpha,E=E),function(x){return(sum(abs(ExponValues-(x[1]+x[2]*x[3]^ExponCycles))))})\$par > > > nls function: > > Emodel<-try(nls(ExponValues ~ (Y0 + a*(E^ExponCycles)) , start= Esp, > algorithm="port"),silent=TRUE) > -- > View this message in context: http://n4.nabble.com/Starting-estimates-for-nls-Exponential-Fit-tp932230p932230.html> Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > [hidden email] mailing list > 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. > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 In reply to this post by Anto Kate Mullen showed one approach to this problem by using DEOptim to get some good starting values. However, I believe that the real issue is scaling (Warning: well-ridden hobby-horse!). With appropriate scaling, as below, nls does fine. This isn't saying nls is perfect -- I've been trying to figure out how to do a nice job of helping folk to scale their problems. Ultimately, it would be nice to has an nls version that will do the scaling and also watch for some other situations that give trouble. Cheers, JN ## JN test rm(list=ls()) ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,                  2468.32,2778.47) ExponCycles <- c(17,18,19,20,21,22,23,24,25) mod <- function(x) x[1] + x[2]*x[3]^ExponCycles modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) fun <- function(x) sum((ExponValues-mod(x))^2) pa<-c(1,2,3) lo<-c(0,0,0) up<-c(20,20,20) names(pa) <- c("Y0", "a", "E") ## fit w/port and GN Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)),              start=pa, algorithm='port', lower=lo, upper=up,              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa,              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l") ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 You used starting values:    pa <- c(1,2,3) but both algorithms (port and Gauss-Newton) fail if you use the slightly different values:    pa <- c(1,2,3.5) Scaling does not fix the underlying sensitivity to starting values. pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both alg. also fail if there is insufficient spread (e.g., both fail for pa <- c(1,1.25,1.5) ). For the record, DE uses (by default at least) a random start and for bad starts will sometimes fail to give good results; decrease the probability this happens by increasing itermax from the default itermax=200, as in: ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),               control=list(trace=FALSE, itermax=1000)) On Wed, 2 Dec 2009, Prof. John C Nash wrote: > Kate Mullen showed one approach to this problem by using DEOptim to get > some good starting values. > > However, I believe that the real issue is scaling (Warning: well-ridden > hobby-horse!). > > With appropriate scaling, as below, nls does fine. This isn't saying nls > is perfect -- I've been trying to figure out how to do a nice job of > helping folk to scale their problems. Ultimately, it would be nice to > has an nls version that will do the scaling and also watch for some > other situations that give trouble. > > Cheers, JN > > > ## JN test > rm(list=ls()) > > ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, >                  2468.32,2778.47) > > ExponCycles <- c(17,18,19,20,21,22,23,24,25) > > mod <- function(x) x[1] + x[2]*x[3]^ExponCycles > > modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) > > fun <- function(x) sum((ExponValues-mod(x))^2) > > > > pa<-c(1,2,3) > lo<-c(0,0,0) > up<-c(20,20,20) > names(pa) <- c("Y0", "a", "E") > > ## fit w/port and GN > Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), >              start=pa, algorithm='port', lower=lo, upper=up, >              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) > > Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, >              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) > > ## fit > matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l") > > ______________________________________________ > [hidden email] mailing list > 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. > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 Kate is correct. The parameter scaling helps quite a bit, but not enough to render the problem "nice" so that many "reasonable" starting points will give useful results. Indeed, a run using "all.methods=TRUE" in our optimx package (on r-forge at http://r-forge.r-project.org/R/?group_id=395) gives                           par  fvalues      method   fns grs itns conv 4   2.194931, 1.000001, 1.000027 566407.6     nlm    NA  NA   3    0 1 2.1949335, -9.0413113, 0.7516435 566407.6   BFGS    37   6 NULL    0 2  2.1950009, 0.2548318, 1.1163498 566404.6 Nelder    82  NA NULL    0 3     1.974226, 1.829957, 1.681338 2409.771   SANN 10000  NA NULL    0 6  1.9904045, 0.6151421, 1.7559401 1696.497  ucminf    51  51 NULL    0 5  1.9906488, 0.6012996, 1.7575365 1696.248  nlminb    80 166   51    0    KKT1  KKT2 4 FALSE FALSE 1  TRUE FALSE 2 FALSE FALSE 3 FALSE  TRUE 6 FALSE  TRUE 5 FALSE  TRUE > A bit of a dog's dinner. On the positive side, this is a simple but nasty problem to illustrate lots of the issues with nonlinear parameter estimation. JN Katharine Mullen wrote: > You used starting values: >    pa <- c(1,2,3) > > but both algorithms (port and Gauss-Newton) fail if you use the slightly > different values: >    pa <- c(1,2,3.5) > > Scaling does not fix the underlying sensitivity to starting values. > pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both > alg. also fail if there is insufficient spread (e.g., both fail for > pa <- c(1,1.25,1.5) ). > > For the record, DE uses (by default at least) a random start and for bad > starts will sometimes fail to give good results; decrease the probability > this happens by increasing itermax from the default itermax=200, as in: > > ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10), >               control=list(trace=FALSE, itermax=1000)) > > On Wed, 2 Dec 2009, Prof. John C Nash wrote: > >> Kate Mullen showed one approach to this problem by using DEOptim to get >> some good starting values. >> >> However, I believe that the real issue is scaling (Warning: well-ridden >> hobby-horse!). >> >> With appropriate scaling, as below, nls does fine. This isn't saying nls >> is perfect -- I've been trying to figure out how to do a nice job of >> helping folk to scale their problems. Ultimately, it would be nice to >> has an nls version that will do the scaling and also watch for some >> other situations that give trouble. >> >> Cheers, JN >> >> >> ## JN test >> rm(list=ls()) >> >> ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, >>                  2468.32,2778.47) >> >> ExponCycles <- c(17,18,19,20,21,22,23,24,25) >> >> mod <- function(x) x[1] + x[2]*x[3]^ExponCycles >> >> modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) >> >> fun <- function(x) sum((ExponValues-mod(x))^2) >> >> >> >> pa<-c(1,2,3) >> lo<-c(0,0,0) >> up<-c(20,20,20) >> names(pa) <- c("Y0", "a", "E") >> >> ## fit w/port and GN >> Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), >>              start=pa, algorithm='port', lower=lo, upper=up, >>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) >> >> Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, >>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) >> >> ## fit >> matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l") >> >> ______________________________________________ >> [hidden email] mailing list >> 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. >> ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 Figuring out the best parameterization for this kind of model is a bit tricky until you get the hang of it. Let the function be         y_t = y_0 + alpha * E^t where uppercase Y_t denotes an observed value and lower case y_t is a predicted value. Index the times by t_1 .... t_n WLOG assume that t_1 is the smallest time and t_n is the largest time. Now we already have a good idea what the predicted values  y_1 and y_n should be as we have observations for them. We have the two equations        y_1 = y_0 + alpha * E^t_1        y_n = y_0 + alpha * E^t_n we can solve these for alpha in terms of y_1,y_n,and E giving      alpha = (y_n-y_1)/(E^t_n -E^t_1)       (1) and solve for y_0  as          y_0 = y_1 - alpha * E^t_1  using (1) for alpha Now use the good estimates Y_1 and Y_n  as the starting values for y_1 and y_n and try some "reasonable value for E (say 0.1 < E < 100) Do the minimization in two stages first holding y_1 and y_n fixed and then estimate y_1,y_n,and E together. This converges in less than a second using AD Model Builder and the starting values (large value for E. 2018.34 2778.47 exp(10.0)  where I deliberately used exp(10) as a large initial value for E. The parameter estimates together with the est std devs are   1   y_1     1.9994e+03 3.9177e-01   2   y_n     2.7881e+03 6.7557e-01   3   log_E  5.6391e-01 1.2876e-03   4   alpha      6.0130e-04 1.9398e-05   5   y_0     1.9906e+03 4.5907e-01   6   E      1.7575e+00 4.3935e-02 There are problems for E near 1 which need to be dealt with if you have to go there, but that is just a technicality. This idea also works well for a logistic say 3 4 or 5 parameter. -- David A. Fournier P.O. Box 2040, Sidney, B.C. V8l 3S3 Canada Phone/FAX 250-655-3364 http://otter-rsch.com______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 In reply to this post by Prof J C Nash (U30A) Thanks everybody, This has been quite helpful, the problem remains tricky but at least now I've got a version of my script that handles all my reactions without error. The DEoptim solution produced good starting values for a lot of reactions, but sadly not for all. I now use scaled parameters and allow more iterations than the standard 50 (1000). I doubt it is completely stable and I think I will run into reactions that will fail the fit, but for now everything works fine so I won't be complaining. Thanks a lot, Antoon   Prof. John C Nash wrote Kate is correct. The parameter scaling helps quite a bit, but not enough to render the problem "nice" so that many "reasonable" starting points will give useful results. Indeed, a run using "all.methods=TRUE" in our optimx package (on r-forge at http://r-forge.r-project.org/R/?group_id=395) gives                           par  fvalues      method   fns grs itns conv 4   2.194931, 1.000001, 1.000027 566407.6     nlm    NA  NA   3    0 1 2.1949335, -9.0413113, 0.7516435 566407.6   BFGS    37   6 NULL    0 2  2.1950009, 0.2548318, 1.1163498 566404.6 Nelder    82  NA NULL    0 3     1.974226, 1.829957, 1.681338 2409.771   SANN 10000  NA NULL    0 6  1.9904045, 0.6151421, 1.7559401 1696.497  ucminf    51  51 NULL    0 5  1.9906488, 0.6012996, 1.7575365 1696.248  nlminb    80 166   51    0    KKT1  KKT2 4 FALSE FALSE 1  TRUE FALSE 2 FALSE FALSE 3 FALSE  TRUE 6 FALSE  TRUE 5 FALSE  TRUE > A bit of a dog's dinner. On the positive side, this is a simple but nasty problem to illustrate lots of the issues with nonlinear parameter estimation. JN Katharine Mullen wrote: > You used starting values: >    pa <- c(1,2,3) > > but both algorithms (port and Gauss-Newton) fail if you use the slightly > different values: >    pa <- c(1,2,3.5) > > Scaling does not fix the underlying sensitivity to starting values. > pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both > alg. also fail if there is insufficient spread (e.g., both fail for > pa <- c(1,1.25,1.5) ). > > For the record, DE uses (by default at least) a random start and for bad > starts will sometimes fail to give good results; decrease the probability > this happens by increasing itermax from the default itermax=200, as in: > > ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10), >               control=list(trace=FALSE, itermax=1000)) > > On Wed, 2 Dec 2009, Prof. John C Nash wrote: > >> Kate Mullen showed one approach to this problem by using DEOptim to get >> some good starting values. >> >> However, I believe that the real issue is scaling (Warning: well-ridden >> hobby-horse!). >> >> With appropriate scaling, as below, nls does fine. This isn't saying nls >> is perfect -- I've been trying to figure out how to do a nice job of >> helping folk to scale their problems. Ultimately, it would be nice to >> has an nls version that will do the scaling and also watch for some >> other situations that give trouble. >> >> Cheers, JN >> >> >> ## JN test >> rm(list=ls()) >> >> ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, >>                  2468.32,2778.47) >> >> ExponCycles <- c(17,18,19,20,21,22,23,24,25) >> >> mod <- function(x) x[1] + x[2]*x[3]^ExponCycles >> >> modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) >> >> fun <- function(x) sum((ExponValues-mod(x))^2) >> >> >> >> pa<-c(1,2,3) >> lo<-c(0,0,0) >> up<-c(20,20,20) >> names(pa) <- c("Y0", "a", "E") >> >> ## fit w/port and GN >> Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), >>              start=pa, algorithm='port', lower=lo, upper=up, >>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) >> >> Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, >>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) >> >> ## fit >> matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l") >> >> ______________________________________________ >> R-help@r-project.org mailing list >> 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 https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 I thought maybe my suggestion for reparameterizing this simple problem was ignored because I didn't supply R code for the problem. Here it is using optim for the optimization. It converges trivially with an initial value for E of 1000. As I stated before, there is nothing at all difficult about this problem. You simply need to parameterize it properly. Of course that is not to say that rescaling is not useful as well, but the important thing is to parameterize the model properly. ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27 78.47) Expon=c(17,18,19,20,21,22,23,24,25) # Example starting estimate calculation E=1000.0 y1=2018 yn=2778.47 nobs=9 #keep y1 and yn fixed and get initial value for E Esp1 <- optim(c(E=E),method ="BFGS",   function(x)   {     E=x[1]     a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])     Y0=y1-a*E^Expon[1];     diff=ExponValues-(Y0+a*E^ExponCycles)     return(1000*sum(diff*diff))   })\$par E=Esp1[1] Esp <- optim(c(y1=y1,yn=yn,E=E),method ="BFGS",   function(x)   {     E=x[3]     a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])     Y0=x[1]-a*E^Expon[1];     diff=ExponValues-(Y0+a*E^ExponCycles)     return(1000*sum(diff*diff))   })\$par y1=Esp[1] y2=Esp[2] E=Esp[3] a=(y2-y1)/(E^Expon[nobs]-E^Expon[1]) Y0=y1-a*E^Expon[1]; -- David A. Fournier P.O. Box 2040, Sidney, B.C. V8l 3S3 Canada Phone/FAX 250-655-3364 http://otter-rsch.com______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Starting estimates for nls Exponential Fit

 Thanks to Dennis Murphy who pointed out that ExponCycles is undefined. It is an R gotcha. I had shortened the name but R still remembered it so the script worked but only on my computer. This should fix that. ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27 78.47) Expon=c(17,18,19,20,21,22,23,24,25) # Example starting estimate calculation E=1000.0 y1=2018 yn=2778.47 nobs=9 #keep y1 and yn fixed and get initial value for E Esp1 <- optim(c(E=E),method ="BFGS",   function(x)   {     E=x[1]     a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])     Y0=y1-a*E^Expon[1];     diff=ExponValues-(Y0+a*E^Expon)     return(1000*sum(diff*diff))   })\$par E=Esp1[1] Esp <- optim(c(y1=y1,yn=yn,E=E),method ="BFGS",   function(x)   {     E=x[3]     a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])     Y0=x[1]-a*E^Expon[1];     diff=ExponValues-(Y0+a*E^Expon)     return(1000*sum(diff*diff))   })\$par y1=Esp[1] y2=Esp[2] E=Esp[3] a=(y2-y1)/(E^Expon[nobs]-E^Expon[1]) Y0=y1-a*E^Expon[1]; -- David A. Fournier P.O. Box 2040, Sidney, B.C. V8l 3S3 Canada Phone/FAX 250-655-3364 http://otter-rsch.com______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.