This below is not solvable with uniroot to find "a": fn=function(a){ b=(0.7/a)-a (1/(a+b+1))-0.0025 } uniroot(fn,c(-500,500)) gives "Error in uniroot(fn, c(-500, 500)) : f() values at end points not of opposite sign" I read R-help posts and someone wrote a function: http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92407.html but it is not very precise. Is there any '"standard" function in R that can solve this? thanks. -- View this message in context: http://www.nabble.com/How-to-solve-difficult-equations--tf3643595.html#a10175603 Sent from the R help mailing list archive at Nabble.com.
I don't see the problem, except that you might want to think about what the error message is telling you. A little exploration of your function always helps, too. > ss <- seq(-2, 2, len = 100) > plot(ss, fn(ss), type = "l") > uniroot(fn, c(-1, 1)) Erreur dans uniroot(fn, c(-1, 1)) : f() values at end points not of opposite sign > fn(-1) [1] 3.330833 > fn(0) [1] -0.0025 > fn(1) [1] 0.5857353 > uniroot(fn, c(-1, 0)) $root [1] -0.6999466 $f.root [1] -13118.83 $iter [1] 18 $estim.prec [1] 7.70751e-05 > uniroot(fn, c(0, 1)) $root [1] 0.001760625 $f.root [1] 8.86832e-06 $iter [1] 3 $estim.prec [1] 6.103516e-05> This below is not solvable with uniroot to find "a": > fn=function(a){ > b=(0.7/a)-a > (1/(a+b+1))-0.0025 > } > uniroot(fn,c(-500,500)) gives > "Error in uniroot(fn, c(-500, 500)) : f() values at end points not of > opposite sign" > > I read R-help posts and someone wrote a function: > http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92407.html > but it is not very precise. Is there any '"standard" function in R > that can > solve this? thanks.
plot(fn,-1,1) could be helpful, hth, Ingmar On 25 Apr 2007, at 09:15, francogrex wrote:> > This below is not solvable with uniroot to find "a": > fn=function(a){ > b=(0.7/a)-a > (1/(a+b+1))-0.0025 > } > uniroot(fn,c(-500,500)) gives > "Error in uniroot(fn, c(-500, 500)) : f() values at end points not of > opposite sign" > > I read R-help posts and someone wrote a function: > http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92407.html > but it is not very precise. Is there any '"standard" function in R > that can > solve this? thanks. > -- > View this message in context: http://www.nabble.com/How-to-solve- > difficult-equations--tf3643595.html#a10175603 > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > 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.Ingmar Visser Department of Psychology, University of Amsterdam Roetersstraat 15 1018 WB Amsterdam The Netherlands t: +31-20-5256735 [[alternative HTML version deleted]]
On 25-Apr-07 07:15:55, francogrex wrote:> > This below is not solvable with uniroot to find "a": > fn=function(a){ > b=(0.7/a)-a > (1/(a+b+1))-0.0025 > } > uniroot(fn,c(-500,500)) gives > "Error in uniroot(fn, c(-500, 500)) : f() values at end points > not of opposite sign" > > I read R-help posts and someone wrote a function: > http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92407.html > but it is not very precise. Is there any '"standard" function > in R that can solve this? thanks.Two answers: Yes, and No. First, "No": Let alpha denote 0.0025, and beta 0.7 (in your function "fn"). Then fna <- function(alpha,beta){ beta*alpha/(1 - alpha) } solves it. But this is not a standard R function. Second, "Yes": and the standard R function is uniroot(). But you can only apply it usefully if you first study the behaviour of your function fn(), in rather careful detail. Over your range (-500,500): a<-10*(-50:50) plot(a,fn(a),pch="+") Clearly something extreme happens just to the left of a=0. So: a <- 0.025*(-100:0) plot(a,fn(a),pch="+") and so for this set of values of 'a' the previous behaviour cannot be seen. So: a <- 0.01*(-100:100)+0.001 plot(a,fn(a),pch="+") so the function goes very negative somewhere around a = -0.7. But fn(500) [1] 0.996102 so it is positive for a=500. Now find (inspired by the latest plot): a[which(fn(a) < (-100))] [1] -0.699 and now you can use uniroot: uniroot(fn,c(-0.699,500)) $root [1] 0.001771128 $f.root [1] 2.379763e-05 $iter [1] 16 $estim.prec [1] 6.103516e-05 and, if that doesn't look precise enough: uniroot(fn,c(-0.699,500),tol=1e-10) $root [1] 0.001754386 $f.root [1] 1.354602e-14 $iter [1] 18 $estim.prec [1] 5e-11 Now compare with the function fna() that solves it directly: fna(0.0025,0.7) [1] 0.001754386 (so in fact it was worth increasing the precision for uniroot). But the lesson to be drawn from all this is that for functions like fn(), which have singularities (here at a = -0.7), the blind application of root-finding functions may not work, since they are not set up to explore the function is the kind of way illustrated above. While there are procedures in the numerical analysis world to handle this kind of thing, they tend to be written for particular classes of function, and again you will have to do a bit of exploration to find out which function to use. And (while someone more knowledgeable may well disagree with me) I suspect that these are not "standard" R funnctions. Hoping this is helpful, Ted. -------------------------------------------------------------------- E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861 Date: 25-Apr-07 Time: 09:31:29 ------------------------------ XFMail ------------------------------
At 03:15 AM 4/25/2007, francogrex wrote:>This below is not solvable with uniroot to find "a": >fn=function(a){ >b=(0.7/a)-a >(1/(a+b+1))-0.0025 >} >uniroot(fn,c(-500,500)) gives >"Error in uniroot(fn, c(-500, 500)) : f() values at end points not of >opposite sign" > >I read R-help posts and someone wrote a function: >http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92407.html >but it is not very precise. Is there any '"standard" function in R that can >solve this? thanks.Actually, if you're solving fn(a)==0, then some trivial algebra leads to a linear equation with a=0.001754. Why are you trying to solve this numerically? Is it a single instance of a larger, more general problem? ===============================================================Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: ral at lcfltd.com Least Cost Formulations, Ltd. URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239 Fax: 757-467-2947 "Vere scire est per causas scire"