Your script is missing something (in particular kw). I presume you are trying to estimate the pK values. You may have more success with package nlsr than nls(). nlsr::nlxb() tries to get the Jacobian of the model specified by a formula and do so by applying symbolic or automatic differentiation. The multi expression function would probably not work. (That is one advantage of nls(), but it has an unstabilized solver and unless carefully called will use a simple forward derivative approximation.) There is also nlsr::nlfb() that can use functions for the residuals and jacobian. Your function is messy, but likely the jacobian can be developed with a bit of work.
Some of the symbolic derivative features of R can help here. Alternatively, there are several numerical approximations in nlsr and the vignette "Intro to nlsr" explains how to use these. Note that simple forward and backward approximations are generally not worth using, but central approximation is quite good. The solver in the nlsr package allows of bounds on the parameters. Since you want them ordered, you may want to transform pK1 < pK2 <- pK3 to pk1, deltapk2, deltapk3 so pK2 = pk1+deltapk2 and pk3 = pk2 + deltapk3 = pk1 + deltapk2 + deltapk3 and all values bounded below by 0 or possibly some small numbers to keep the paramters apart. I won't pretend any of this is trivial. It is VERY easy to make small errors that still allow for output that is disastrously incorrect. If it is possible to get a single "formula" as one expression even if spread over multiple lines, then nlxb() might be able to handle it. J Nash (maintainer of nlsr and optimx) On 2023-11-06 11:53, Troels Ring wrote: ....
HEPESFUNC <- function(H,SID,HEPTOT,pK1,pK2,pK3) { XX <- (H^3/(10^-pK3*10^-pK2*10^-pK1)+H^2/(10^-pK3*10^-pK2)+H/(10^-pK3)) IV <- HEPTOT/(1+XX) I <- IV*H^3/(10^-pK3*10^-pK2*10^-pK1) II <- IV*H^2/(10^-pK3*10^-pK2) III <- IV*H/10^-pK3 H - kw/H + SID + I*2 + II - IV } HEPTOT <- 0.050 SID <- c(-seq(10,1)*1e-4,0,seq(1,10)*1e-4) pHobs <- c(4.63,4.68,4.72,4.77,4.83,4.9,4.96,5.04,5.12,5.21,5.3, 5.39,5.48,5.55,5.63,5.69,5.74,5.8,5.85,5.89,5.93) pK1 <- -1; pK2 <- 3; pK3 <- 7.55 # literature values pK3 <- 7.6; pK2 <- 2.96 # values eye-balled to be better pH <- c() for (i in 1:length(SID)) { pH[i] <- -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4, HEPTOT=HEPTOT,SID = SID[i], pK1=pK1,pK2=pK2,pK3=pK3)$root) } plot(SID,pH) points(SID,pHobs,col="red")
______________________________________________ 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.