[R] R and Excel disagreement - Goal Seek versus uniroot
Peter Dalgaard
p.dalgaard at biostat.ku.dk
Sat Apr 12 11:13:12 CEST 2008
Troels Ring wrote:
> Dear friends - occurring in Windows R2.6.2
> I am modeling physical chemistry in collaboration with a friend who has
> preferred working in Excel. I used uniroot, and find a solution to a two
> buffer problem in acid-base chemistry which I believe is physiologically
> sensible. Using "goal seek" in Excel my friend found another plausible
> root, quite close to zero, and a plot of the function fd below shows
> that it may be about OK - but Excel fails to find the plausible root
> indicated by uniroot.
>
> fd <- function(H,SID,KA1,KA2,ATOT1,ATOT2){
> H^4+H^3*(SID+KA1+KA2)+H^2*(SID*(KA1+KA2)+KA1*KA2-ATOT1*KA1-ATOT2*KA2-kw)+
> H*(SID*KA1*KA2-KA1*KA2*(ATOT1+ATOT2)-kw*(KA1+KA2))-kw*KA1*KA2}
>
> KA1 <- 1e-6;KA2 <- 1e-5;ATOT1 <- 0.05; ATOT2 <- 0.03; SID <-
> 0.05;kw <- 4.4e-14
> options(digits=20)
> uniroot(fd,c(0,1),tol =
> .Machine$double.eps,maxiter=100000,KA1=KA1,KA2=KA2,SID=SID,
> ATOT1=ATOT1,ATOT2=ATOT2)$root #1.16219184891115e-06
>
>
> And here is the plot indicating to me that the root I found at 1.16e-06
> might be right and also that the root at 0 might be OK too - but less
> interesting.
>
> H <- seq(0,2e-6,length=10000)
> plot(H,fd(H,SID,KA1,KA2,ATOT1,ATOT2))
> abline(v=1.16219184891115e-06)
> abline (h=0)
> text(9.0e-7,3e-19,"H=1.1622e-6")
>
>
> However, using optimize another solution is found:
>
> xmin <- optimize(fd,c(0,1),tol =
> .Machine$double.eps,KA1=KA1,KA2=KA2,SID=SID,
> ATOT1=ATOT1,ATOT2=ATOT2)
> xmin
> #$minimum
> #[1] 6.102746508205e-07
> #$objective
> #[1] -9.72253673562293e-20
>
>
> I would very much appreciate some help on doing this modeling with very
> small numbers in R. Comments on the capability of Excel's goal seek
> would also be welcome.
>
> Best wishes
> Troels
>
>
A number of points to be made here:
(1) uniroot finds a root, optimize a minimum. In this case, your
function is fairly close to a quadratic so the plot looks like a
parabola with two root and a minimum about midways. (And a good starting
point could be to drop off the 3rd and 4th order terms in H and analyze
the quadratic using your highschool arithmetic.)
(2) Since the value at zero is -4.4e-25 and the function is decreasing
there, the root on the left is not zero, but slightly to the left of it
(at -1.47e-12, it would seem), so the only root in the interval (0,1) is
indeed the one at 1.16e-6. Presumably, Excel is using local
linearization around zero and finding the nearest root.
(3) you _can_ use optimize (or optim) to find roots, but you need to
square the objective function first, and you need rather better starting
values, because the squared function is not convex (plus you have the
risk of finding an optimum which is not at zero).
(4) you have to be very careful with optimizers and root finders if the
scale of the objective function and its parameters is not close to 1,
because they can (and will) misinterpret the small values as indications
of convergence. For uniroot and optimize, you may need to set the
tolerance well below .Machine$double.eps (this is floating point, so
double.eps is only relevant for values that are not themselves small),
and for optim. you can play with control parameters fnscale and
parscale, e.g.
> uniroot(fd,c(-1e-6,0),tol =
+ 1e-70,maxiter=100000,KA1=KA1,KA2=KA2,SID=SID,
+ ATOT1=ATOT1,ATOT2=ATOT2)
$root
[1] -1.466662866313071e-12
$f.root
[1] 0
$iter
[1] 4
$estim.prec
[1] 3.187027620317954e-19
> optimize(fq,c(-1e-6,0),tol =
+ 1e-70,KA1=KA1,KA2=KA2,SID=SID,
+ ATOT1=ATOT1,ATOT2=ATOT2)
$minimum
[1] -1.466662866319826e-12
$objective
[1] 4.10609522514202e-72
> optim(0, fq, control=list(fnscale=1e-70, parscale=1e-12),method="BFGS",
+ KA1=KA1,KA2=KA2,SID=SID,
+ ATOT1=ATOT1,ATOT2=ATOT2)
$par
[1] -1.466662866312404e-12
$value
[1] 4.000708456650843e-74
$counts
function gradient
510 20
$convergence
[1] 0
$message
NULL
(5) This is a fourth degree polynomial in H. R has a polyroot() function...
--
O__ ---- Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
More information about the R-help
mailing list