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

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

G'day Troels, On Mon, 6 Nov 2023 20:43:10 +0100 Troels Ring <tring at gvdnet.dk> wrote:> Thanks a lot! This was amazing. I'm not sure I see how the conditiion > pK1 < pK2 < pK3 is enforced?One way of enforcing such constraints (well, in finite computer arithemtic only "<=" can be enforced) is to rewrite the parameters as: pK1 = exp(theta1) ## only if pK1 > 0 pK2 = pK1 + exp(theta2) pK3 = pk2 + exp(theta3) And then use your optimiser to optimise over theta1, theta2 and theta3. There might be better approaches depending on the specific problem. Cheers, Berwin