0

I receive an error from nls function in R. I search some similar questions, but do not solve this problem. For example, I try to use nlsLM from library 'minpack.lm', it also fails. So I have to ask for help here. Following is the code:

tt = c(10, 30, 50, 90, 180, 360, 720, 1440, 2880, 4320, 8640, 12960)
x = c(
  1.53901e-06,
  1.22765e-06,
  1.11200e-06,
  9.25185e-07,
  8.71809e-07,
  8.80705e-07,
  8.36225e-07,
  7.82849e-07,
  8.18433e-07,
  6.04928e-07,
  3.46944e-07,
  4.44800e-07
)
y = c(
  3.81639e-06,
  5.00623e-06,
  4.62815e-06,
  5.10631e-06,
  4.48359e-06,
  3.30487e-06,
  2.64879e-06,
  2.13727e-06,
  8.02865e-07,
  1.91487e-06,
  3.73855e-06,
  2.32631e-06
)

nt = length(tt)
L0 = 0.005
y0 = 0.000267681

model = function(K, Kd, k1) {
  
  eta = 5 / (4 * Kd + 40)
  eta1 = 1 - eta
  eta1_seq = eta1 ^ c(0:(nt - 1))
  Lt = L0 * eta * cumsum(eta1_seq)
  
  b = K * x - K * Lt + 1
  L = (-b + sqrt(b ^ 2.0 + 4 * K * Lt)) / (2 * K)
  
  cx = x * K * L / (K * L + 1)
  qx = Kd * cx
  
  q1 = y0 * (1 - k1 * sqrt(tt))
  
  y = qx + q1
  
  return(y)
  
}

fit <- nls(
  y ~ model(K, Kd, k1),
  start = list(K = 1e+15,
               Kd = 10,
               k1 = 1e-5),
  lower = c(1e+13, 1, 1e-10),
  upper = c(1e+20, 200, 1e-3),
  algorithm = "port"
)

Thanks in advance for your help!

Lin
  • 195
  • 1
  • 12
  • You need better starting values (will only work if your data supports the model at all). – Roland Oct 30 '20 at 07:55
  • It is hard to figure out the starting values. I sample 3000 rows of values, and choose the one with minimum error as the starting values, but I still receive same error. – Lin Oct 30 '20 at 08:11
  • There is an `alpha` in your `model` and in `start` that does not do anything. It has 0 gradient everywhere, so an error is thrown. If you remove it and set `K = 1e3` as a starting value, at least the model fits. – Bas Oct 30 '20 at 08:44
  • The starting value of K is too low to the real problem. I edit the code with lower and upper limits. – Lin Oct 30 '20 at 09:03
  • Are you sure in the `sqrt` there shouldn't be `- 4 * K * Lt` instead of `+`? – Bas Oct 30 '20 at 09:22
  • Yes, in the sqrt function it is "b ^ 2.0 + 4 * K * Lt"! – Lin Oct 30 '20 at 13:04
  • 1
    try `plot(y ~ tt); lines(yy ~ tt, col = "red")` for a variety of fitted values parameter combinations where yy is the fitted values from a particular set of parameters. If you can't find any parameter sets whose fitted values yy look reasonable then there may be no value of the parameters that represent your data. Also virtually all nonlinear optimization algorithms require that the parameters be of similar order of magnitude so you will need to scale your parameters to achieve this if in the first step you were able to find some reasonable combinations. – G. Grothendieck Oct 30 '20 at 15:02

0 Answers0