Search code examples
rregressioncurve-fittingnls

Exponential regression using nls


I am having trouble fitting an exponential curve to my data.

Here's my code:

x<-c(0.134,0.215,0.345,0.482,0.538,0.555)
y<-c(0,0,0.004,0.291,1.135,1.684)
plot(x,y)
estimates<- list(b1 = 0.1, b2 = 5e-7)
nlfit <- nls(y ~ b1 * (exp(x/b2)-1), start=estimates)
lines(x, predict(nlfit), col = 2)

But I get the following error:

Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model

I have tried several approaches described on Stack Overflow, such as removing the zeros or fitting it to a simpler model (log), but all of them just gave me a different error. My guess would be that I need better starting values, but I can't seem to find them without getting any errors.


Solution

  • Just guessing list(b1=1, b2=1) works fine. (Filling in all ones is a reasonable default/desperation strategy if all of the predictor variables in your model are reasonably scaled ...)

    nlfit <- nls(y ~ b1 * (exp(x/b2)-1), start=list(b1=1,b2=1))
    lines(x, predict(nlfit), col = 2)
    

    It's a little bit hard to do the usual trick of converting to a log-linear model in order to find good starting estimates because you have the -1 term modifying the exponential ... However, you can try to figure out a little more about the geometry of the curve and eyeball the data you have:

    • at x=0, y should be approximately equal to 0 (OK, that doesn't give us much in the way of information about parameter values)
    • b2 represents the "e-folding time" of the curve, which looks like it's on the order of about 0.1 or maybe a little bit less (the curve looks like it increases about 5x between x=0.4 and x=0.5, but that's about the right order of magnitude)
    • the y-value at x=0.5 is about 0.5, so a reasonable starting value for b1 is approximately 0.5/(exp(0.5/0.1)-1) = 0.003.

    The final estimates are b1=3.1e-6, b2=4.2e-2, but the starting values were close enough to give a reasonably sensible fit, which is usually good enough.

    nlfit <- nls(y ~ b1 * (exp(x/b2)-1), start=list(b1=1,b2=1), data=data.frame(x,y))
    plot(x,y)
    curve(0.003*(exp(x/0.1)-1), col="blue", add=TRUE)
    xvec <- seq(0.1,0.6,length=51)
    lines(xvec,predict(nlfit,newdata=data.frame(x=xvec)),col="red")
    

    points, initial guess,and final fit