#StackBounty: #r #non-linear-regression Find initial conditions for nonlinear models using the nlsLM function

Bounty: 50

I am using the nlsLM function from the minpack.lm package to find the values of parameters a, e, and c that give the best fit to the data out.
Here is my code:

n <- seq(0, 70000, by = 1)
TR <- 0.946
b <- 2000
k <- 50000
nr <- 25
na <- 4000
nd <- 3200
d <- 0.05775

y <- d + ((TR*b)/k)*(nr/(na + nd + nr))*n
## summary(y)
out <- data.frame(n = n, y = y)
plot(out$n, out$y)

## Estimate the parameters of a nonlinear model
library(minpack.lm)
k1 <- 50000
k2 <- 5000

fit_r <- nlsLM(y ~ a*(e*n + k1*k2 + c), data=out,
               start=list(a = 2e-10,
                          e = 6e+05, 
                          c = 1e+07), lower = c(0, 0, 0), algorithm="port")
print(fit_r)
## summary(fit_r)

df_fit <- data.frame(n = seq(0, 70000, by = 1))
df_fit$y <- predict(fit_r, newdata = df_fit)
plot(out$n, out$y, type = "l", col = "red", ylim = c(0,10))
lines(df_fit$n, df_fit$y, col="green")
legend(0,ceiling(max(out$y)),legend=c("observed","predicted"), col=c("red","green"), lty=c(1,1), ncol=1)

The fitting to the data seems to be very sensitive to initial conditions. For example:

  • with list(a = 2e-10, e = 6e+05, c = 1e+07), this gives a good fit:
Nonlinear regression model
  model: y ~ a * (e * n + k1 * k2 + c)
   data: out
        a         e         c 
2.221e-10 5.895e+05 9.996e+06 
 residual sum-of-squares: 3.225e-26

Algorithm "port", convergence message: Relative error between `par' and the solution is at most `ptol'.

enter image description here

  • with list(a = 2e-01, e = 100, c = 2), this gives a bad fit:
Nonlinear regression model
  model: y ~ a * (e * n + k1 * k2 + c)
   data: out
        a         e         c 
1.839e-08 1.000e+02 0.000e+00 
 residual sum-of-squares: 476410

Algorithm "port", convergence message: Relative error in the sum of squares is at most `ftol'.

enter image description here

So, Is there an efficient way to find initial conditions that give a good fit to the data ?

EDIT:

I added some conditions to define initial conditions for the parameters a, e, and c:

Using the result of the linear model lm(y ~ n):

Call:
lm(formula = y ~ n, data = out)

Coefficients:
(Intercept)            n  
  0.0577500    0.0001309
c = 0.0577500/a - k1*k2 > 0 and 
e = 0.0001309/a > 0
0 < a < 0.0577500/(k1*k2)

With this, is it possible to find the values of parameters that give the best fit to the data ?


Get this bounty!!!

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.