r-bisect

# Differences

This shows you the differences between two versions of the page.

 — r-bisect [2017/09/20 19:05] (current) 2012/05/06 21:29 matt created 2012/05/06 21:29 matt created Line 1: Line 1: + ======Bisection root-finding in R====== + =====Iterative===== + The following code was adapted from an original version by Ravi Varadhan posted on the R-help mailing list 20/9/2010. + + <​code>​ + 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) + } + ​ + + =====Recursive===== + The following code was posted to R-help by myself (BioStatMatt) as a recursive alternative to the iterative method above. + + <​code>​ + 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, ...)) + } + ​ + + =====Illustration===== + + These functions find the root of a smooth function, provided two '​bracketing'​ starting values: + + <​code>​ + > testFn <- function(x, a = 1) exp(-x) - a*x + > curve(testFn,​ from = -1, to = 3, xlab = "​x",​ ylab = "​testFn(x)"​) + > # Plot the '​bracketing points'​ + > points(c(0, 2),​testFn(c(0,​2)),​ pch = "​x",​ col= "​red"​ ) + > # Compute and plot the root + > abline(h = 0, lty = 2, col = "​red"​) + > abline(v = bisectMatt(testFn,​ 0, 2), lty = 2, col = "​red"​) + ​ + + This code produces the following graphic: + + {{:​bisect.png|Bisection}} + + =====Timing===== + The iterative function is roughly twice faster than the recursive function: + <​code>​ + > testFn <- function(x, a) exp(-x) - a*x + > system.time(replicate(1000,​ bisectMatt(testFn,​ 0, 2, a=1))) + ​user ​ system elapsed ​ + 0.656   ​0.000 ​  ​0.658 ​ + > system.time(replicate(1000,​ bisectRavi(testFn,​ 0, 2, a=1))) + ​user ​ system elapsed ​ + 0.384   ​0.000 ​  ​0.393 ​ +
r-bisect.txt ยท Last modified: 2017/09/20 19:05 (external edit)