Search code examples
rnls

R: nls2 misses the solution

I'm trying to fit a bi exponential function:

t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))
ggplot(df, aes(t, y)) + geom_line() +  scale_y_continuous(limits=c(0, 50))

enter image description here

This problem can't be solve by a simple transformation like log so I wanted to use the nls2 package:

library(nls2)

fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)
fit <- nls2(fo,
            start = list(Ahat=5, Bhat=5, alphahat=0.5,betahat=0.5),
            algorithm = "brute-force",
            trace = TRUE,
            lower = c(Ahat=0, Bhat=0, alphahat=0, betahat=0),
            upper = c(Ahat=50, Bhat=50, alphahat=10,betahat=10))
fit

Here is the result:

Nonlinear regression model
  model: y ~ Ahat * exp(-alphahat * t) + Bhat * exp(-betahat * t)
   data: NULL
    Ahat     Bhat alphahat  betahat 
     5.0      5.0      0.5      0.5 
 residual sum-of-squares: 37910

Number of iterations to convergence: 4 
Achieved convergence tolerance: NA

I assume something is wrong in my code because :

  • data: NULL ?
  • Why only 4 iterations ?
  • Hard to think nls2 didn't find a better solution than the starting point.
  • The result is far from the solution
Solution

  • From the documentation, the start parameter should be a data.frame of two rows that define the grid to search in, or a data.frame with more rows corresponding to parameter combinations to test if you are using brute-force. Also, nls will have trouble with your fit because it is a perfect curve, there is no noise. The brute-force method is slow, so here is an example where the search space is decreased for nls2. The result of the brute-force nls2 is then used as the starting values with nls default algorithm (or you could use nls2), after adding a tiny bit of noise to the data.

    ## Data
t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))

library(nls2)
fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)

## Define the grid to search in,
## Note: decreased the grid size
grd <- data.frame(Ahat=c(10,30),
                  Bhat=c(10, 30),
                  alphahat=c(0,2),
                  betahat=c(0,1))

## Do the brute-force
fit <- nls2(fo,
            data=df,
            start = grd,
            algorithm = "brute-force",
            control=list(maxiter=100))
coef(fit)
#       Ahat       Bhat   alphahat    betahat 
# 10.0000000 23.3333333  0.0000000  0.3333333 

## Now, run through nls:
## Fails, because there is no noise
final <- nls(fo, data=df, start=as.list(coef(fit)))

## Add a little bit of noise
df$y <- df$y+rnorm(nrow(df),0,0.001)
coef((final <- nls(fo, data=df, start=as.list(coef(fit)))))
#        Ahat        Bhat    alphahat     betahat 
# 10.00034000 19.99956016  0.01000137  0.25000966 

## Plot
plot(df, col="steelblue", pch=16)
points(df$t, predict(final), col="salmon", type="l")

    enter image description here