Trying to match data with R and nls by function with condition in it
I am trying to fit some data to a function that has a validity constraint in it. More precisely, a function with a different value if t <= T and t> T.
Here is the code I 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 am getting the following error:
Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
Is it because my initial parameters are not good enough (I've tried several), my problem is poorly translated to R, or is it a fundamental math mistake I missed?
source to share
nls(...)
uses the default Gaussian Newton method; this error message, which is quite common in fact, means that the jacobian matrix cannot be inverted.
I think your problem is related to the fact that your composite function (RHS of your formula) is not continuous in t=tswitch
for arbitrary values of other parameters. To put it another way, the requirement for the function to be continuous puts a constraint on other parameters - they are independent of each other. Also, the derivative of a compound function will never be continuous at t=tswitch
- yours posExpDecay(...)
has a positive derivative for all t
, while yours negExpDecay(...)
has a negative derivative for all t
.
I cannot know if there is a theoretical reason for this functional form, but these +/- exponents are usually modeled using the product of positive and negative decay as shown below.
Note. I usually use nlsLM(...)
a package minpack.lm
that uses the much more robust Levenberg – Marquardt algorithm. It has the same signature as the function nls(...)
in the R base.
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 well enough.
source to share