[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)
-log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4,
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
minpack.lm:
minpack.lm::nlsLM(
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,
Ivan
More information about the R-help
mailing list