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.

Reply via email to