I was just reading about the merge sort algorithm last night (BTW, here
is a fun link http://www.youtube.com/watch?v=t8g-iYGHpEA). There are
some interesting similarities in this context. Here's a recursive method
for bisection:
bisectMatt <- function(fn, lo, hi, tol = 1e-7, ...) {
flo <- fn(lo, ...)
fhi <- fn(hi, ...)
if(flo * fhi > 0)
stop("root is not bracketed by lo and hi")
mid <- (lo + hi) / 2
fmid <- fn(mid, ...)
if(abs(fmid) <= tol || abs(hi-lo) <= tol)
return(mid)
if(fmid * fhi > 0)
return(bisectMatt(fn, lo, mid, tol, ...))
return(bisectMatt(fn, mid, hi, tol, ...))
}
# Adapted from Ravi's original
bisectRavi <- function(fn, lo, hi, tol = 1e-7, ...) {
flo <- fn(lo, ...)
fhi <- fn(hi, ...)
if (flo * fhi > 0)
stop("root is not bracketed by lo and hi")
chg <- hi - lo
while (abs(chg) > tol) {
mid <- (lo + hi) / 2
fmid <- fn(mid, ...)
if (abs(fmid) <= tol) break
if (flo * fmid < 0) hi <- mid
if (fhi * fmid < 0) lo <- mid
chg <- hi - lo
}
return(mid)
}
testFn <- function(x, a) exp(-x) - a*x
> system.time(bM <- bisectMatt(testFn, 0, 2, a=1))
user system elapsed
0.000 0.000 0.001
> system.time(bR <- bisectRavi(testFn, 0, 2, a=1))
user system elapsed
0.000 0.000 0.001
> bM
[1] 0.5671433
> bR
[1] 0.5671433
Of course, Ravi's version is better for production (and most likely
faster, though not significantly so in this example) because recursion
is more expensive than looping.
-Matt
On Fri, 2010-09-17 at 17:44 -0400, Ravi Varadhan wrote:
> Here is something simple (does not have any checks for bad input), yet
> should be adequate:
>
> bisect <- function(fn, lower, upper, tol=1.e-07, ...) {
> f.lo <- fn(lower, ...)
> f.hi <- fn(upper, ...)
> feval <- 2
>
> if (f.lo * f.hi > 0) stop("Root is not bracketed in the specified interval
> \n")
> chg <- upper - lower
>
> while (abs(chg) > tol) {
> x.new <- (lower + upper) / 2
> f.new <- fn(x.new, ...)
> if (abs(f.new) <= tol) break
> if (f.lo * f.new < 0) upper <- x.new
> if (f.hi * f.new < 0) lower <- x.new
> chg <- upper - lower
> feval <- feval + 1
> }
> list(x = x.new, value = f.new, fevals=feval)
> }
>
> # An example
> fn1 <- function(x, a) {
> exp(-x) - a*x
> }
>
> bisect(fn1, 0, 2, a=1)
>
> bisect(fn1, 0, 2, a=2)
>
>
> Ravi.
>
> -----Original Message-----
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf Of Peter Dalgaard
> Sent: Friday, September 17, 2010 4:16 PM
> To: Gregory Gentlemen
> Cc: r-help@r-project.org
> Subject: Re: [R] Is there a bisection method in R?
>
> On 09/17/2010 09:28 PM, Gregory Gentlemen wrote:
> > If uniroot is not a bisection method, then what function in R does use
> bisection?
> >
>
> Why do you assume that there is one? uniroot contains a better algorithm
> for finding bracketed roots.
>
> It shouldn't be too hard to roll your own if you need one for
> pedagogical purposes.
>
--
Matthew S. Shotwell
Graduate Student
Division of Biostatistics and Epidemiology
Medical University of South Carolina
______________________________________________
R-help@r-project.org 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.
