R uniroot 一維求根(零)

說明

函數uniroot在從lower到upper的區間中搜索函數f相對於其第一個參數的根(即零)。

將extendInt設置為非"no"字符串，意味著如果sign(f(x))在區間終點不滿足要求，則尋找正確的interval = c(lower,upper)；請參閱“詳細信息”部分。

用法

uniroot(f, interval, ...,
        lower = min(interval), upper = max(interval),
        f.lower = f(lower, ...), f.upper = f(upper, ...),
        extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,
        tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)

參數

f	求根的函數。
interval	包含要搜索根的區間的 end-points 的向量。
...	要傳遞給 f 的其他命名或未命名參數
lower , upper	要搜索的區間的下端點和上端點。
f.lower , f.upper	分別與 f(upper) 和 f(lower) 相同。一旦 f() 包含重要的計算，從通常已知的調用者處傳遞這些值會更經濟。
extendInt	字符串，指定當 f() 在端點處沒有不同的符號時，是否應擴展間隔 c(lower,upper) 或直接產生錯誤。默認值 "no" 保留搜索間隔，因此會產生錯誤。可以縮寫。
check.conv	邏輯指示底層 uniroot 的收斂警告是否應被捕獲為錯誤，以及 maxiter 迭代中的不收斂是否應為錯誤而不是警告。
tol	所需的精度(收斂容差)。
maxiter	最大迭代次數。
trace	整數；如果是肯定的，則產生追蹤信息。更高的值提供更多細節。

細節

請注意，... 之後的參數必須完全匹配。

必須指定interval 或同時指定lower 和upper：上端點必須嚴格大於下端點。對於 extendInt="no" 來說，端點處的函數值必須具有相反的符號(或零)，這是默認值。否則，如果 extendInt="yes" ，則在兩側擴展間隔，以搜索符號更改，即，直到搜索間隔 [l,u] 滿足 f(l) \cdot f(u) \le 0 。

如果知道f如何在根x_0處改變符號，也就是說，如果函數在那裏增加或減少，則extendInt可以(通常應該)指定為"upX"(對於“upward crossing”)或"downX" ，分別。同樣，定義 S := \pm 1 ，以在解決方案中要求 S = \mathrm{sign}(f(x_0 + \epsilon)) 。在這種情況下，搜索間隔[l,u]可能被擴展為使得S\cdot f(l)\le 0和S \cdot f(u) \ge 0。

uniroot() 使用基於下麵參考文獻中給出的算法的 Fortran 子例程 zeroin(來自 Netlib)。它們假設一個連續函數(已知該函數在區間內至少有一個根)。

如果 f(x) == 0 或算法一步的 x 變化小於 tol(加上 x 中的表示錯誤容限)，則聲明收斂。

如果算法在 maxiter 步驟中未收斂，則會打印警告並返回當前近似值。

f將被稱為f(x, ...)對於數值x.

傳遞給 f 的參數具有特殊語義，用於在調用之間共享。該函數不應該複製它。

值

至少包含五個組件的列表：root 和 f.root 給出根的位置以及在該點計算的函數的值。 iter 和 estim.prec 給出使用的迭代次數以及 root 的近似估計精度。 (如果根出現在端點之一，則估計精度為 NA 。) init.it 包含初始 extendInt 迭代次數(如果有)，否則為 NA 。在這種extendInt迭代的情況下，iter包含其中的sum和zeroin迭代。

將來可能會添加更多組件。

例子

require(utils) # for str

## some platforms hit zero exactly on the first step:
## if so the estimated precision is 2/3.
f <- function (x, a) x - a
str(xmin <- uniroot(f, c(0, 1), tol = 0.0001, a = 1/3))

## handheld calculator example: fixed point of cos(.):
uniroot(function(x) cos(x) - x, lower = -pi, upper = pi, tol = 1e-9)$root

str(uniroot(function(x) x*(x^2-1) + .5, lower = -2, upper = 2,
            tol = 0.0001))
str(uniroot(function(x) x*(x^2-1) + .5, lower = -2, upper = 2,
            tol = 1e-10))

## Find the smallest value x for which exp(x) > 0 (numerically):
r <- uniroot(function(x) 1e80*exp(x) - 1e-300, c(-1000, 0), tol = 1e-15)
str(r, digits.d = 15) # around -745, depending on the platform.

exp(r$root)     # = 0, but not for r$root * 0.999...
minexp <- r$root * (1 - 10*.Machine$double.eps)
exp(minexp)     # typically denormalized


##--- uniroot() with new interval extension + checking features: --------------

f1 <- function(x) (121 - x^2)/(x^2+1)
f2 <- function(x) exp(-x)*(x - 12)

try(uniroot(f1, c(0,10)))
try(uniroot(f2, c(0, 2)))
##--> error: f() .. end points not of opposite sign

## where as  'extendInt="yes"'  simply first enlarges the search interval:
u1 <- uniroot(f1, c(0,10),extendInt="yes", trace=1)
u2 <- uniroot(f2, c(0,2), extendInt="yes", trace=2)
stopifnot(all.equal(u1$root, 11, tolerance = 1e-5),
          all.equal(u2$root, 12, tolerance = 6e-6))

## The *danger* of interval extension:
## No way to find a zero of a positive function, but
## numerically, f(-|M|) becomes zero :
u3 <- uniroot(exp, c(0,2), extendInt="yes", trace=TRUE)

## Nonsense example (must give an error):
tools::assertCondition( uniroot(function(x) 1, 0:1, extendInt="yes"),
                       "error", verbose=TRUE)

## Convergence checking :
sinc <- function(x) ifelse(x == 0, 1, sin(x)/x)
curve(sinc, -6,18); abline(h=0,v=0, lty=3, col=adjustcolor("gray", 0.8))

uniroot(sinc, c(0,5), extendInt="yes", maxiter=4) #-> "just" a warning


## now with  check.conv=TRUE, must signal a convergence error :

uniroot(sinc, c(0,5), extendInt="yes", maxiter=4, check.conv=TRUE)


### Weibull cumulative hazard (example origin, Ravi Varadhan):
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u

n <- 1000
u <- -log(runif(n))
a <- 1/2
b <- 1
## Find failure times
ru <- sapply(u, function(x)
   uniroot(froot, u=x, a=a, b=b, interval= c(1.e-14, 1e04),
           extendInt="yes")$root)
ru2 <- sapply(u, function(x)
   uniroot(froot, u=x, a=a, b=b, interval= c(0.01,  10),
           extendInt="yes")$root)
stopifnot(all.equal(ru, ru2, tolerance = 6e-6))

r1 <- uniroot(froot, u= 0.99, a=a, b=b, interval= c(0.01, 10),
             extendInt="up")
stopifnot(all.equal(0.99, cumhaz(r1$root, a=a, b=b)))

## An error if 'extendInt' assumes "wrong zero-crossing direction":

uniroot(froot, u= 0.99, a=a, b=b, interval= c(0.1, 10), extendInt="down")

來源

基於 'zeroin.c' 在https://netlib.org/c/brent.shar.

參考

Brent, R. (1973) Algorithms for Minimization without Derivatives. Englewood Cliffs, NJ: Prentice-Hall.

也可以看看

polyroot 用於多項式的所有複數根； optimize、nlm。