Skip to content

NLS fit for exponential distribution

5 messages · Diviya Smith, Peter Dalgaard, Achim Zeileis +1 more

Original
Peter Dalgaard
#
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.

            

      
      
      
    

        
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.

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Achim Zeileis
    

    

      
      #
    

  


  

      
On Sun, 12 Jun 2011, peter dalgaard wrote:

        

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      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: