? 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. -- Best regards, Ivan
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. >