[R] non-linear regression and root finding

[J C Nash]

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 <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.
>>
> 
> ______________________________________________
> R-help using 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.