Search code examples
rnonlinear-optimizationnls

nls() : “Error singular gradient matrix at initial parameter estimates ”

I've read many similar questions but still couldn't find the answer. Here is some data that I'm using to calibrate the equation below:

set.seed(100)
i <- sort(rexp(n = 100,rate = 0.01))
Tr <- sort(runif(n = 100,min = 5,max = 100))

k_start <- 3259
u_start <- 0.464
t0_start <- 38
n_start <- -1

i_test <- k_start*Tr^u_start * (5 + t0_start)^n_start

m <- nls(i~(k * Tr^u / (5+t0)^n), start = list(k = k_start, u = u_start,
                                               t0 = t0_start, n = n_start))

When I used nlsLM and the same error came up:

Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates

For the start values, I tried to use the values from calibration in Python and still the same error occurs.

There's also another way to use that equation that is like this: However, the result is the same error.

d_start <- 43

m <- nls(i ~ (k * Tr^u / d),
         start = list(k = k_start, u = u_start,d=d_start))

When I use only the numerator it works, but that's not what I need. Any help will be very much appreciated.

Solution

  • In the first nls, the right hand side depends on k, t0 and n only through k / (5+t0)^n so it is over parameterized as one parameter could represent their combined effect. In the second nls the right hand side depends only on k and d through k / d so again the problem has been over parameterized and one parameter could represent their combined effect.

    Getting rid of the excess parameters and getting the starting values using a linear model it converges.

    fit.lm <- lm(log(i) ~ log(Tr))
co <- coef(fit.lm)
fit <- nls(i ~ k * Tr ^ u, start = list(k = exp(co[[1]]), u = co[[2]]))
fit
## Nonlinear regression model
##   model: i ~ k * Tr^u
##    data: parent.frame()
##         k         u 
## 0.0002139 3.0941602 
##  residual sum-of-squares: 79402
##
## Number of iterations to convergence: 43 
## Achieved convergence tolerance: 5.354e-06

    Reciprocal Model

    Below we fit a "reciprocal model" which has the same number of parameters but a better fit as measured by the deviance which is the residual sum of squares. A lower value means better fit.

    # reciprocal model
fit.recip <- nls(i ~ 1/(a + b * log(Tr)), start = list(a = 1, b = 1))

deviance(fit)
## [1] 79402.17
deviance(fit.recip)
## [1] 25488.1

    Graphics

    Below we plot both fit (red) and fit.recip (blue) models.

    plot(i ~ Tr)
lines(fitted(fit) ~ Tr, col = "red")
lines(fitted(fit.recip) ~ Tr, col = "blue")
legend("topleft", legend = c("fit", "fit.recip"), lty = 1, col = c("red", "blue"))

    (continued after plot) screenshot

    plinear

    Note that the plinear algorithm could be used as an alternative algorithm to fit the fit model above to avoid having to supply a starting value for k. It also has the additional benefit that it requires substantially fewer iterations in this case (14 vs. 45). With plinear the formula should omit the linear argument, k, as it is implied by the algorithm and will be reported as .lin .

    nls(i ~ Tr ^ u, start = list(u = co[[2]]), algorithm = "plinear")
## Nonlinear regression model
##   model: i ~ Tr^u
##    data: parent.frame()
##         u      .lin 
## 3.0941725 0.0002139 
##  residual sum-of-squares: 79402
##
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 3.848e-06