使用 uniroot 查找根目录

Question

我正在尝试使用 uniroot() 函数查找以下函数（基于 Gamma (gamma()) 函数）的根：

cv = 0.056924/1.024987^2

fx2 = function(theta, eta){
  p1 = 1 - 2/(theta*(1-eta))
  p2 = 1 - 1/(theta*(1-eta))
  return(( gamma(p1)/(gamma(p2))^2 ) - (cv+1) )
}

这个函数给出了以下情节：

v = seq(0, 1, 0.01)
plot(v, fx2(3.0, v), type='l' )

我觉得这个函数的根接近0.33，但是uniroot()函数没有找到根，返回结果如下：

uniroot(fx2, interval = c(0,0.3), theta=3 )

Error in uniroot(fx2, interval = c(0, 0.3), theta = 3) : f() values at end points not of opposite sign

如何找到这个函数的根？有没有其他包有更准确的算法？

Answer 1

我首先将您的函数重写为（可选）根据首先在对数刻度上（通过 lgamma()）然后求幂的计算来表达 gamma(p1)/gamma(p2)^2。这在数值上更稳定，结果将在下面变得清楚......（我可能搞砸了对数尺度计算 - 你应该仔细检查它。Update/warning：更仔细地阅读文档 (!!)，lgamma() 求值为 gamma 函数的 绝对值 的对数。所以可能会出现一些奇怪的符号在下面的答案中。事实仍然是，如果您正在评估 x<0 的伽马函数的比率（即在该值可能变为负值的情况下），很可能会发生 Bad Stuff。

cv = 0.056924/1.024987^2
fx3 <- function(theta, eta, lgamma = FALSE) {
  p1 <- 1 - 2/(theta*(1-eta))
  p2 <- 1 - 1/(theta*(1-eta))
  if (lgamma) {
    val <- exp(lgamma(p1) - 2*lgamma(p2)) - (cv+1)
  } else {
    val <- ( gamma(p1)/(gamma(p2))^2 ) - (cv+1)
  }
}

使用和不使用对数缩放计算函数：

x <- seq(0, 1, length.out = 20001)
v <- sapply(x, fx3, theta = 3.0, lgamma =  TRUE)
v2 <- sapply(x, fx3, theta = 3.0, lgamma =  FALSE)

查找根（下面有更多解释）：

uu <- uniroot(function(eta) fx3(3.0, eta, lgamma = TRUE),
        c(0.4, 0.5))

绘制它：

par(las=1, bty="l")
plot(x, abs(v), col = as.numeric(v<0) + 1, type="p", log="y",
     pch=".", cex=3)
abline(v = uu$root, lty=2)
cvec <- sapply(c("blue","magenta"), adjustcolor, alpha.f = 0.2)
points(x, abs(v2), col=cvec[as.numeric(v2<0) + 1], pch=".", cex=3)

这里我在对数刻度上绘制了 绝对值 ，符号用颜色表示（black/blue >0，红色>洋红色 <0）。 Black/red是log-scale计算，blue/magenta是原始计算。我还以非常高的分辨率绘制了该函数，以尽量避免丢失或错误描述特征。

这里发生了很多奇怪的事情。

两个版本的函数在 x=1/3 附近做了一些有趣的事情；原始版本看起来像一个极点（值从-∞发散到+∞，“returns”），而对数尺度计算上升到+∞和returns而不改变符号。
对数尺度计算在 x=0.45 附近有一个根（符号翻转时绝对值变小），但原始计算没有——大概是因为某种灾难性的精度损失？如果我们给 uniroot 个不包括极点的边界，它可以找到这个根。
在更大的 x 值处还有更多的极点 and/or 根，我没有探索过。

所有这些基本上都表明，在不知道其数学属性是什么的情况下乱用这个函数是非常危险的。我通过数值探索发现了一些东西，但最好分析函数，这样你才能真正知道发生了什么；如果函数的行为足够奇怪，任何数值探索都可能被愚弄。

使用 uniroot 查找根目录

Find root using uniroot

r

function

uniroot