I won't send to list, but just to the two of you, as I don't have
anything to add at this time. However, I'm wondering if this approach
is worth writing up, at least as a vignette or blog post. It does need
a shorter example and some explanation of the "why" and some testing
perhaps.
If there's interest, I'll be happy to join in. And my own posting suggests
how the ordering is enforced by bounding the "delta" parameters from below.
Note that nls() can only handle bounds in the "port" algorithm, and the
man page rather pours cold water on using that algorithm.
Best, JN
On 2023-11-06 14:43, Troels Ring wrote:
Thanks a lot! This was amazing. I'm not sure I see how the conditiion pK1 < pK2 < pK3 is enforced? - it comes from the
derivation via generalized Henderson-Hasselbalch but perhaps it is not really necessary. Anyway, the use of Vectorize
did the trick!
Best wishes
Troels
Den 06-11-2023 kl. 19:19 skrev Ivan Krylov:
В Mon, 6 Nov 2023 17:53:49 +0100
Troels Ring <tr...@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.
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
