[R] non-linear regression and root finding

Ivan Krylov kry|ov@r00t @end|ng |rom gm@||@com
Mon Nov 6 19:19:02 CET 2023

В Mon, 6 Nov 2023 17:53:49 +0100
Troels Ring <tring using gvdnet.dk> пишет:

> Hence I wonder if I could somehow have non linear regression to find
> the 3 pK values. Below is HEPESFUNC which delivers charge in the
> fluid for known pKs, HEPTOT and SID. Is it possible to have
> root-finding in the formula with nls?

Sure. Just reformulate the problem in terms of a function that takes a
vector of predictors (your independent variable SID) and the desired
parameters (pK1, pK2, pK3) as separate arguments and then returns
predicted values of the dependent variable (to compare against pHobs):

kw <- 1e-14 # I'm assuming
pHm <- Vectorize(function(SID, pK1, pK2, pK3)
        HEPTOT=HEPTOT,SID = SID, pK1=pK1,pK2=pK2,pK3=pK3)$root))

(Yes, Vectorize() doesn't make the function any faster, but I don't see
an easy way to rewrite this function to make it truly vectorised.)

Unfortunately, nls() seems to take a step somewhere where crossprod()
of the Jacobian of the model function cannot be inverted and fails with
the message "Singular gradient". I wish that R could have a more
reliable built-in nonlinear least squares solver. (I could also be
holding it wrong.) Meanwhile, we have excellent CRAN packages nlsr and

 pHobs ~ pHm(SID, pK1, pK2, pK3),
 data.frame(pHobs = pHobs, SID = SID),
 start = c(pK1 = pK1, pK2 = pK2, pK3 = pK3),
 # the following is also needed to avoid MINPACK failing to fit
 lower = rep(-1, 3), upper = rep(9, 3)
# Nonlinear regression model
#   model: pHobs ~ pHm(SID, pK1, pK2, pK3)
#    data: data.frame(pHobs = pHobs, SID = SID)
#    pK1     pK2    pK3
# -1.000   2.966  7.606
#  residual sum-of-squares: 0.001195
# Number of iterations to convergence: 15 
# Achieved convergence tolerance: 1.49e-08

(Unfortunately, your code seemed to have a lot of non-breakable spaces,
which confuse R's parser and make it harder to copy&paste it.)

I think that your function can also be presented as a degree-5
polynomial in H, so it should also be possible to use polyroot() to
obtain your solutions in a more exact manner.

Best regards,

More information about the R-help mailing list