1

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

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' )

曲线图显示在 x ~ 0.3 处发散的值

在我看来，这个函数的根接近 0.33，但是该uniroot()函数没有找到根，返回以下结果：

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

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

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

4

1 回答 1

3

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

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，红色>洋红色<0）。黑色/红色是对数尺度计算，蓝色/品红色是原始计算。我还以非常高分辨率绘制了函数，以避免丢失或错误描述特征。

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

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

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

于 2021-06-17T21:43:37.327 回答