r-bisect
Differences
This shows you the differences between two versions of the page.
|
r-bisect [2017/09/20 19:05] |
r-bisect [2017/09/20 19:05] (current) |
||
|---|---|---|---|
| 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) | ||
| + | } | ||
| + | </code> | ||
| + | |||
| + | =====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, ...)) | ||
| + | } | ||
| + | </code> | ||
| + | |||
| + | =====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") | ||
| + | </code> | ||
| + | |||
| + | 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 | ||
| + | </code> | ||
r-bisect.txt ยท Last modified: 2017/09/20 19:05 (external edit)
Page Tools
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 4.0 International
