[R] non-linear regression and root finding

[Troels Ring]

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 <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.
>