[R] NLS fit for exponential distribution

[Dennis Murphy]

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 <diviya.smith at gmail.com> 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]]
>
> ______________________________________________
> 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.
>