R › R help

NLS fit for exponential distribution

Classic

List

Threaded

5 messages Options

Diviya Smith

NLS fit for exponential distribution

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

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 . Any suggestion how I can do this... Any help would be most
appreciated. Thanks in advance.

Diviya

[[alternative HTML version deleted]]

______________________________________________
[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.

Peter Dalgaard-2

Re: NLS fit for exponential distribution

On Jun 12, 2011, at 18:57 , Diviya Smith wrote:

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? 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: [hidden email] Priv: [hidden email]

______________________________________________
[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.

Achim Zeileis-4

Re: NLS fit for exponential distribution

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: [hidden email] Priv: [hidden email]
>
> ______________________________________________
> [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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

djmuseR

Re: NLS fit for exponential distribution

In reply to this post by Diviya Smith

Hi:

If you use RStudio, then you can use its manipulate package to figure
out starting values for the model visually through sliders. Walmes
Zeviani posted the template of how to do this last week; I've just
adapted it for the models under consideration. It's interesting to
play with this because it shows how strongly some of the parameters
are correlated with one another...and it's fun :) Thanks to Walmes for
the excellent tutorial example.

Firstly, x and y (or the inputs in general) need to be in the global
environment as vectors of the same length.
[RStudio is a GUI, so I'm assuming that you run things from its
command line.] Load the manipulate package (which comes with RStudio)
and use the manipulate() function to create a plot of the data and fit
a curve to it. The sliders adjust the parameter values and you play
with them until the curve fits closely to the data. The range of the
sliders should match the range of plausible values you think they
might take. The output of the manipulate() function is a vector of
starting values that you can pass into nls(), which is done by closing
the window containing the sliders.

# Model 1:

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)
plot(x, y) # Initial plot of the data
start <- list() # Initialize an empty list for the starting values

library(manipulate)
manipulate(
{
plot(x, y)
k <- kk; b0 <- b00; b1 <- b10
curve(k*exp(-b1*x) + b0, add=TRUE)
start <<- list(k=k, b0=b0, b1=b1)
},
kk=slider(0.01, 0.2, step = 0.01, initial = 0.1),
b10=slider(0.01, 0.4, step = 0.01, initial = 0.3),
b00=slider(-0.01, -0.001, step=0.001,initial= -0.005))

# When done, start() is a list of named parameters in
# the global environment
# Model fit:
fit1 <- nls(y ~ k*exp(-b1*x) + b0, start = start)
summary(fit1)

### Model 2: [following Peter Dalgaard's suggested model]
### Use the estimates from fit1 and shrink b1 to anticipate
### potential effect of b2; the initial estimate in the slider for
### b1 will be too big, so it needs to be shrunk (a lot :)
start <- list()
manipulate(
{
plot(x, y)
k <- kk; b0 <- b00; b1 <- b10; b2 <- b20
curve(k*(exp(-b1*x) + exp(-b2*x)) + b0, add=TRUE)
start <<- list(k=k, b0=b0, b1=b1, b2 = b2)
},
kk=slider(0.01, 0.1, step = 0.01, initial = 0.04),
b10=slider(0.01, 0.4, step = 0.01, initial = 0.3),
b20 = slider(0.01, 0.1, step = 0.005, initial = 0.05),
b00=slider(-0.01, -0.001, step=0.001,initial= -0.004))

fit2 <- nls(y ~ k*(exp(-b1*x) + exp(-b2*x)) + b0, start = start)
summary(fit2)

anova(fit1, fit2)

HTH,
Dennis

On Sun, Jun 12, 2011 at 9:57 AM, Diviya Smith <[hidden email]> 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
>
> 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 . Any suggestion how I can do this... Any help would be most
> appreciated. Thanks in advance.
>
> Diviya
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [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.
>

Diviya Smith

Re: NLS fit for exponential distribution

Thank you for very informative answers. This is great help. Peter, thanks
for correcting the model. That was exactly what I meant - apologies for the
typo.

I was able to run the analysis on some simulated data and the subset of data
that I had posted earlier. However, when I apply the analysis to full
dataset, I usually get the following errors-

fm2 <- nls(y ~ b + a * (exp(m1 * -x) + exp(m2 * -x)), start = list(a = 1, b
= 0, m1 = 0.1, m2 = 0.5), trace = TRUE)
84.42045 : 1.0 0.0 0.1 0.5
0.1593364 : 0.072399953 -0.000226868 0.101135423 0.571436093
0.09032837 : 7.897930e-02 7.905359e-06 1.271435e-01 1.872206e+00
0.06834495 : 0.0808842841 0.0005815865 0.1725302473 4.9048478439
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model

OR

Error in nls(y ~ b + a * (exp(m1 * -x) + exp(m2 * -x)), start = list(a = 1,
:
number of iterations exceeded maximum of 50

Is there a way to optimize the start values? I was unable to try the RStudio
solution as RStudio keeps crashing on my computer.

Also is there any way to set the start value in NLS - I mean if I want to
fit the model to only a subset of data starting at x = 10, how can I specify
that?

Thanks,
Diviya

On Sun, Jun 12, 2011 at 7:19 PM, Dennis Murphy <[hidden email]> wrote:

> Hi:
>
> If you use RStudio, then you can use its manipulate package to figure
> out starting values for the model visually through sliders. Walmes
> Zeviani posted the template of how to do this last week; I've just
> adapted it for the models under consideration. It's interesting to
> play with this because it shows how strongly some of the parameters
> are correlated with one another...and it's fun :) Thanks to Walmes for
> the excellent tutorial example.
>
> Firstly, x and y (or the inputs in general) need to be in the global
> environment as vectors of the same length.
> [RStudio is a GUI, so I'm assuming that you run things from its
> command line.] Load the manipulate package (which comes with RStudio)
> and use the manipulate() function to create a plot of the data and fit
> a curve to it. The sliders adjust the parameter values and you play
> with them until the curve fits closely to the data. The range of the
> sliders should match the range of plausible values you think they
> might take. The output of the manipulate() function is a vector of
> starting values that you can pass into nls(), which is done by closing
> the window containing the sliders.
>
> # Model 1:
>
> 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)
> plot(x, y) # Initial plot of the data
> start <- list() # Initialize an empty list for the starting values
>
> library(manipulate)
> manipulate(
> {
> plot(x, y)
> k <- kk; b0 <- b00; b1 <- b10
> curve(k*exp(-b1*x) + b0, add=TRUE)
> start <<- list(k=k, b0=b0, b1=b1)
> },
> kk=slider(0.01, 0.2, step = 0.01, initial = 0.1),
> b10=slider(0.01, 0.4, step = 0.01, initial = 0.3),
> b00=slider(-0.01, -0.001, step=0.001,initial= -0.005))
>
> # When done, start() is a list of named parameters in
> # the global environment
> # Model fit:
> fit1 <- nls(y ~ k*exp(-b1*x) + b0, start = start)
> summary(fit1)
>
> ### Model 2: [following Peter Dalgaard's suggested model]
> ### Use the estimates from fit1 and shrink b1 to anticipate
> ### potential effect of b2; the initial estimate in the slider for
> ### b1 will be too big, so it needs to be shrunk (a lot :)
> start <- list()
> manipulate(
> {
> plot(x, y)
> k <- kk; b0 <- b00; b1 <- b10; b2 <- b20
> curve(k*(exp(-b1*x) + exp(-b2*x)) + b0, add=TRUE)
> start <<- list(k=k, b0=b0, b1=b1, b2 = b2)
> },
> kk=slider(0.01, 0.1, step = 0.01, initial = 0.04),
> b10=slider(0.01, 0.4, step = 0.01, initial = 0.3),
> b20 = slider(0.01, 0.1, step = 0.005, initial = 0.05),
> b00=slider(-0.01, -0.001, step=0.001,initial= -0.004))
>
> fit2 <- nls(y ~ k*(exp(-b1*x) + exp(-b2*x)) + b0, start = start)
> summary(fit2)
>
> anova(fit1, fit2)
>
>
> HTH,
> Dennis
>
> On Sun, Jun 12, 2011 at 9:57 AM, Diviya Smith <[hidden email]>
> 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
> >
> > 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 . Any suggestion how I can do this... Any help would be most
> > appreciated. Thanks in advance.
> >
> > Diviya
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > [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.
> >
>

[[alternative HTML version deleted]]

______________________________________________
[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.