swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: 113add4c7a9bfd37466a9b41cf2dc060b7b87091 authored by HwB on 06 December 2012, 00:00:00 UTC
version 1.3.1
version 1.3.1
Tip revision: 113add4
fibsearch.R
##
## f i b s e a r c h . R Fibonacci Search
##
fibsearch <- function(f, a, b, ...,
endp = FALSE, tol = .Machine$double.eps^(1/2))
# Fibonacci search for a univariate function minimum in a bounded interval
{
if (a >= b) stop("Left endpoint a must be smaller than b.")
tol <- max(tol, .Machine$double.eps)
# Compute Fibonacci numbers [F0,] F1, F2, ..., Fm
F <- c(1, 2); n <- 2
while (F[n] <= 2*(b-a)/tol) { F[n+1] <- F[n] + F[n-1]; n <- n + 1 }
# Initialize values (k == 0)
x1 <- a; x2 <- b
xa <- a + (b-a) * F[n-2]/F[n]; fxa <- f(xa, ...)
xb <- a + (b-a) * F[n-1]/F[n]; fxb <- f(xb, ...)
# Compute iteration
k <- 1
while (k <= n-3 && xa < xb && (x2 - x1) >= tol) {
if (fxa > fxb) {
x1 <- xa; xa <- xb
xb <- x1 + (x2-x1) * F[n-k-1]/F[n-k]
fxa <- fxb; fxb <- f(xb, ...)
} else {
x2 <- xb; xb <- xa
xa <- x1 + (x2-x1) * F[n-k-2]/F[n-k]
fxb <- fxa; fxa <- f(xa, ...)
}
k <- k + 1
}
# Finally use the mean and consider endpoints
xmin <- (xa+xb)/2; fmin <- f(xmin, ...)
if (endp) {
fa <- f(a, ...); fb <- f(b, ...)
if (xmin-a < tol && fa < fmin) {
xmin <- a; fmin <- fa
} else {
if (b - xmin < tol && fb < fmin) {
xmin <- b; fmin <- fb
}
}
}
estim.prec <- max(xmin-x1, x2-xmin)
return(list(xmin=xmin, fmin=fmin, niter=k, estim.prec=estim.prec))
}