[R] NLS fit for exponential distribution

[Achim Zeileis]

On Sun, 12 Jun 2011, peter dalgaard wrote:

>
> On Jun 12, 2011, at 18:57 , Diviya Smith wrote:
>
>> Hello there,
>>
>> I am trying to fit an exponential fit using Least squares to some data.
>>
>> #data
>> x <- c(1 ,10,  20,  30,  40,  50,  60,  70,  80,  90, 100)
>> y <- c(0.033823,  0.014779,  0.004698,  0.001584, -0.002017, -0.003436,
>> -0.000006, -0.004626, -0.004626, -0.004626, -0.004626)
>>
>> sub <- data.frame(x,y)
>>
>> #If model is y = a*exp(-x) + b then
>> fit <- nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), trace
>> = TRUE)
>>
>> This works well. However, if I want to fit the model : y = a*exp(-mx)+c then
>> I try -
>> fit <- nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m=
>> 0), trace = TRUE)
>>
>> It fails and I get the following error -
>> Error in nlsModel(formula, mf, start, wts) :
>>  singular gradient matrix at initial parameter estimates
>
>
> If a==0, then a*exp(-m*x) does not depend on m. So don't use a=0 as initial value.
>
>>
>> Any suggestions how I can fix this? Also next I want to try to fit a sum of
>> 2 exponentials to this data. So the new model would be  y = a*exp[(-m1+
>> m2)*x]+c .
>
> That's not a sum of exponentials. Did you mean a*(exp(-m1*x) + 
> exp(-m2*x)) + c?

OK, that makes more sense. Also, scaling of the variables may help. 
Something like this could work:

## scaled data
d <- data.frame(x = c(1, 1:10 * 10), y = 100 *  c(
   0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, -0.000006,
   -0.004626, -0.004626, -0.004626, -0.004626))

## model fits
n1 <- nls(y ~ a*exp(-m * x) + b, data = d,
   start = list(a = 1, b = 0, m = 0.1))
n2 <- nls(y ~ a * (exp(-m1 * x) + exp(-m2 * x)) + b, data = d,
   start = list(a = 1, b = 0, m1 = 0.1, m2 = 0.5))

## visualization
plot(y ~ x, data = d)
lines(d$x, fitted(n1), col = 2)
lines(d$x, fitted(n2), col = 4)

## ANOVA
anova(n1, n2)

## LR test
library("lmtest")
lrtest(n1, n2)

which both seem to indicate that n1 is sufficient.

hth,
Z

> Anyways, same procedure with more parameters. Just 
> beware the fundamental exchangeability of m1 and m2, so don't initialize 
> them to the same value.
>
>> Any suggestion how I can do this... Any help would be most
>> appreciated. Thanks in advance.
>
> -- 
> Peter Dalgaard
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
>
> ______________________________________________
> R-help at 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.
>