r - 使用 uniroot 查找根目录
问题描述
我正在尝试使用该函数找到以下函数的根(基于 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' )
在我看来,这个函数的根接近 0.33,但是该uniroot()
函数没有找到根,返回以下结果:
uniroot(fx2, interval = c(0,0.3), theta=3 )
uniroot(fx2, interval = c(0, 0.3), theta = 3) 中的错误:端点处的 f() 值不是相反的符号
我如何找到这个函数的根?还有其他算法更准确的软件包吗?
解决方案
我首先将您的函数重写为(可选)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 值处还有更多的极点和/或根,我没有探索。
所有这一切基本上都表明,在不知道它的数学属性是什么的情况下乱用这个函数是非常危险的。我通过数值探索发现了一些东西,但最好分析一下函数,这样你才能真正知道发生了什么;如果函数的行为足够奇怪,任何数值探索都可能被愚弄。
推荐阅读
- ros - 脚本执行“休眠”多长时间?
- latex - 在乳胶中使用tabularx时如何将文本包装在单元格内?
- java - 重复时在 XML 文件中找到我想要的标签的最佳方法是什么?
- node.js - 不调用 Jest spiedOn 方法,而是调用一次更改测试用例顺序
- java - 在服务器端和客户端之间共享字段约束
- php - 为什么谷歌地图中多边形的坐标返回错误的纬度、纬度值?
- django - Django ORM:将模型之间的外键关系限制为有限的次数
- javascript - 如何使用 id 和目标 _blank 在新选项卡中打开页面
- python - 如何在 python 中获取 QtLabel 的当前值?- PyQt5
- c++ - 使其成为 C++ 标准的 Boost 库