Nls singular gradient matrix - fitting parameters in integral upper bounds

Question

I am trying to make a binding nlsfor a slightly complex expression that includes two integrals with two matching parameters in their upper limits.

I got an error

"Error in nlsModel (formula, mf, start, wts): singular matrix gradient when evaluating initial parameters."

I already searched in previous answers, but didn't help. Initializing the parameters looks fine, I tried changing the parameters, but no one is working. If my function has only one integral, everything works very well, but when adding the second integral term, an error is simply received. I do not think that the function is overridden, because I performed other tricks with a lot more parameters, and they worked. Below I wrote a list with some data.

Minimal example:

integrand <- function(X) {
  return(X^4/(2*sinh(X/2))^2)
}

fitting = function(T1, T2, N, D, x){
  int1 = integrate(integrand, lower=0, upper = T1)$value
  int2 = integrate(integrand, lower=0, upper = T2)$value
  return(N*(D/x)^2*(exp(D/x)/(1+exp(D/x))^2
)+(448.956*(x/T1)^3*int1)+(299.304*(x/T2)^3*int2))
}
fit = nls(y ~ fitting(T1, T2, N, D, x),
start=list(T1=400,T2=200,N=0.01,D=2))

------> For reference, the work worked:

integrand <- function(X) {
  return(X^4/(2*sinh(X/2))^2)
}
fitting = function(T1, N, D, x){
  int = integrate(integrand, lower=0, upper = T1)$value
  return(N*(D/x)^2*(exp(D/x)/(1+exp(D/x))^2 )+(748.26)*(x/T1)^3*int)
}
fit = nls(y ~ fitting(T1 , N, D, x), start=list(T1=400,N=0.01,D=2))

-------> Data to illustrate the problem:

dat<- read.table(text="x       y
0.38813 0.0198
0.79465 0.02206
1.40744 0.01676
1.81532 0.01538
2.23105 0.01513
2.64864 0.01547
3.05933 0.01706
3.47302 0.01852
3.88791 0.02074
4.26301 0.0256
4.67607 0.03028
5.08172 0.03507
5.48327 0.04283
5.88947 0.05017
6.2988  0.05953
6.7022  0.07185
7.10933 0.08598
7.51924 0.0998
7.92674 0.12022
8.3354  0.1423
8.7384  0.16382
9.14656 0.19114
9.55062 0.22218
9.95591 0.25542", header=TRUE)

I can’t understand what will happen. I need to perform this fit for three integrated components, but even for two I have this problem. I really appreciate your help. Thank.

Answer 1

:

fitting1 <- function(par, x, y) {
  sum((fitting(par[1], par[2], par[3], par[4], x) - y)^2)
}

library(optimx)
res <-  optimx(c(400, 200, 0.01, 2),
       fitting1,
       x = DF$x, y = DF$y,
       control = list(all.methods = TRUE))

print(res)

#                  p1       p2          p3         p4         value fevals gevals niter convcode kkt1 kkt2 xtimes
#BFGS        409.7992 288.6416  -0.7594461   39.00871  1.947484e-03    101    100    NA        1   NA   NA   0.22
#CG          401.1281 210.9087  -0.9026459   20.80900  3.892929e-01    215    101    NA        1   NA   NA   0.25
#Nelder-Mead 414.6402 446.5080  -1.1298606 -227.81280  2.064842e-03     89     NA    NA        0   NA   NA   0.02
#L-BFGS-B    412.4477 333.1338  -0.3650530   37.74779  1.581643e-03     34     34    NA        0   NA   NA   0.06
#nlm         411.8639 333.4776  -0.3652356   37.74855  1.581644e-03     NA     NA    45        0   NA   NA   0.04
#nlminb      411.9678 333.4449  -0.3650271   37.74753  1.581643e-03     50    268    48        0   NA   NA   0.07
#spg         422.0394 300.5336  -0.5776862   38.48655  1.693119e-03   1197     NA   619        0   NA   NA   1.06
#ucminf      412.7390 332.9228  -0.3652029   37.74829  1.581644e-03     45     45    NA        0   NA   NA   0.05
#Rcgmin            NA       NA          NA         NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
#Rvmmin            NA       NA          NA         NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
#newuoa      396.3071 345.1165  -0.3650286   37.74754  1.581643e-03   3877     NA    NA        0   NA   NA   1.02
#bobyqa      410.0392 334.7074  -0.3650289   37.74753  1.581643e-03   7866     NA    NA        0   NA   NA   2.07
#nmkb        569.0139 346.0856 282.6526588 -335.32320  2.064859e-03     75     NA    NA        0   NA   NA   0.01
#hjkb        400.0000 200.0000   0.0100000    2.00000  3.200269e+00      1     NA     0     9999   NA   NA   0.01

Levenberg-Marquardt , nlsLM , nls , :

library(minpack.lm)
fit <- nlsLM(y ~ fitting(T1, T2, N, D, x),
          start=list(T1=412,T2=333,N=-0.36,D=38), data = DF, trace = TRUE)
#It.    0, RSS = 0.00165827, Par. =        412        333      -0.36         38
#It.    1, RSS = 0.00158186, Par. =    417.352    329.978    -0.3652     37.746
#It.    2, RSS = 0.00158164, Par. =    416.397    330.694  -0.365025    37.7475
#It.    3, RSS = 0.00158164, Par. =    416.618    330.568  -0.365027    37.7475
#It.    4, RSS = 0.00158164, Par. =    416.618    330.568  -0.365027    37.7475
#Error in nlsModel(formula, mf, start, wts) : 
#  singular gradient matrix at initial parameter estimates