Search code examples
rregressionnlsnon-linear-regressionmodel-fitting

error fitting function to data using nls


I have some issue using nls() to estimate parameters. I have a following set of functions to explain some data in hand:

funk1 <- function(a,x) { x^2*exp(-(l*(1-exp(-r*a))/r)) }

funk2 <- function(x) { sapply(x, function (s)
{ integrate(funk1, lower = 0, upper = s, x=s)$value }) }

I am trying to fit funk2 to y:

y <- sort(runif(100, 0, 10^8))

When I use nls():

nls(y ~ funk2(z1$days.post.bmt), data= z1, start=list(l=0.02, r=0.002), trace=T)

it shows me following error:

Error in f(x, ...) : object 'l' not found

Isn't the whole point of nls() to substitute different values for parameters l and r from parameter space to fit the function by minimizing SSR and give the parameter estimates? why it needs value of l for it to work? I am definitely missing something big here. Please help!

Thanks in advance!


Solution

  • You must pass parameters l and r as function arguments of funk1 and funk2.

    funk1 <- function(a,x,l,r) {
      x^2*exp(-(l*(1-exp(-r*a))/r))
      }
    
    funk2 <- function(x,l,r) {
      sapply(x, function (s) {
                  integrate(funk1, lower = 0, upper = s, x=s, l=l, r=r)$value
                  })
      }
    

    I will generate some data to test:

    z <- data.frame(days.post.bmt = 1:100,
                    y = funk2(1:100, l = 1, r = 1) + rpois(100, 1:100))
    
    nls(y ~ funk2(days.post.bmt,l,r), data = z, start = list(l = 0.5, r = 0.5))
    
    #Nonlinear regression model
    #  model: y ~ funk2(days.post.bmt, l, r)
    #   data: z
    #     l      r 
    #0.9405 0.9400 
    # residual sum-of-squares: 6709
    
    #Number of iterations to convergence: 5 
    #Achieved convergence tolerance: 2.354e-07
    

    As a counter example, consider:

    bad_funk1 <- function(a,x) {
      x^2*exp(-(l*(1-exp(-r*a))/r))
      }
    
    bad_funk2 <- function(x) {
      sapply(x, function (s) {
                  integrate(funk1, lower = 0, upper = s, x=s)$value
                  })
      }
    
    nls(y ~ bad_funk2(days.post.bmt), data = z, start = list(l = 0.5, r = 0.5))
    # Error in f(x, ...) (from #2) : argument "l" is missing, with no default