使用 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 )

uniroot(fx2, interval = c(0, 0.3), theta = 3) 中的错误： f() 端点的值不是相反的符号

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

Answer 1

我首先将您的函数重写为（可选）表达gamma(p1)/gamma(p2)^2，该计算首先在对数尺度上完成（通过lgamma()）然后取幂。这在数值上更稳定，其后果将在下面变得清晰......（我可能搞砸了对数尺度计算 - 你应该仔细检查它。更新/警告：阅读文档更仔细（！！），lgamma() 评估为 gamma 函数的绝对值的日志。所以下面的答案中可能会出现一些奇怪的迹象。事实仍然是，如果您正在评估 x

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)

这里我在对数刻度上绘制绝对值，符号用颜色表示（黑色/蓝色>0，红色>洋红色

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

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

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