Search code examples
rnls

Trying to fit data with R and nls on a function with a condition in it

I am trying to fit some data to a function which has validity limit in them. More precisely a function with different value if t<=T and t>T.

Here is the code I have tried:

posExpDecay <- function(t,tau,max,toff){ 1+max*(1-exp(-(t-toff)/tau)) } 
negExpDecay <- function(t,tau,max){ 1+max*exp(-(t)/tau) } 

data<-structure(list(t = c(0.67, 1, 1.33, 1.67, 2, 4, 6, 8, 10), y = c(1.02,2.33, 3.08, 3.34, 3.41,2.50, 1.86, 1.44, 1.22)), .Names = c("t", "y"), row.names = c(13L, 17L, 21L, 25L, 29L,37L, 45L, 49L, 53L), class = "data.frame")

fit <- nls(y~ifelse(t<=tswitch,
                    posExpDecay(t,tau1,max1,toff),
                    negExpDecay(t,tau2,max2)),
                  data,
                  start=list(max1=3,tau1=0.7,max2=7,tau2=2,toff=0.1,tswitch=3))

And I get the following error:

Error in nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter estimates

Is this that my starting parameters are not good enough (I tried several), is my problem not well translated in R, or a fundamental mathematical error I missed?

Solution

  • nls(...) uses the Gauss Newton method by default; that error message, which is quite common actually, means that the Jacobian matrix cannot be inverted.

    I think your problem has to do with the fact the your composite function (the RHS of your formula) is not continuous at t=tswitch for arbitrary values of the other parameters. To say it differently, the requirement that the function be continuous puts a constraint on the other parameters - they are not independent of each other. Also, the derivative of the composite function will never be continuous at t=tswitch - your posExpDecay(...) has a positive derivative for all t, whereas your negExpDecay(...) has a negative derivative for all t.

    I can't know if there is a theoretical reason for this functional form, but these +/- exponentials are generally modeled using the product of a positive and negative decay, as shown below.

    Note: I generally use nlsLM(...) in the minpack.lm package, which uses the much more robust Levenberg Marquardt algorithm. It has the same signature as the nls(...) function in base R.

    f <- function(t, max,tau1,tau2,toff) max*exp(-t/tau1)*(1-exp(-(t-toff)/tau2))
library(minpack.lm)
fit <- nlsLM(y~f(t,max,tau1,tau2,toff),data,
             start=list(max=15,tau1=0.7,tau2=2,toff=.2))
summary(fit)
# ...
# Parameters:
#      Estimate Std. Error t value Pr(>|t|)    
# max   4.72907    0.29722  15.911 1.78e-05 ***
# tau1  6.75926    0.54093  12.496 5.82e-05 ***
# tau2  0.51211    0.08209   6.238  0.00155 ** 
# toff  0.53595    0.02667  20.093 5.64e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.113 on 5 degrees of freedom
# 
# Number of iterations to convergence: 19 
# Achieved convergence tolerance: 1.49e-08

plot(y~t,data)
curve(predict(fit,data.frame(t=x)),add=T,col="blue")

    As you can see this much simpler function (fewer parameters) fits reasonably well.