VGAM logit函数溢出

Question

VGAM版本0.93

> logit(1000, inverse=T)
[1] 0 # it should be 1

问题在于：

exp(1000 - log1p(exp(1000)))

在这里log1p(exp(1000))变成了Inf

因此，它使用的数值方法不处理大的数字，相比于plogis在基础上的正确工作。

它是否值得一次错误报告，我可以在哪里提交？

Answer 1

Update：这是一个bug，尽管存在浮点问题，但它应该返回1。实际上，在这种情况下，应该使用plogis函数，因为它正确地处理了问题。作者和维护人员在描述文件中被报道为Thomas，你应该给他发一封电子邮件。

您的机器不能表示那么小的浮点数。考虑逆logit函数：

inv.logit<-function(x) exp(x)/(1+exp(x))

即使是500岁，它也是非常小的：

inv.logit(-500)
# 7.124576e-218

在我自己的机器上，这已经接近机器所能代表的极限了。您可以找到这个值.Machine$double.xmin。

.Machine$double.xmin
# [1] 2.225074e-308 

如果您真的对这里的确切值感兴趣，您将不得不将这些数字转换为可以在您的计算机上表示的比例。

事实上，对于很多人来说，这个问题并没有改变。为了弄清楚您要求机器表示什么(当您要求1000的逆逻辑时)，请尝试使用gmp包。你要求补充这个数字的倒数：

library(gmp)
exp(1)^as.bigz(1000)
Big Integer ('bigz') :
[1] 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376    

问题是这个数字的倒数，它将是非常小的，并且在你的计算机上是无法表示的。

实际上，您可以使用Rmpfr包来计算这个数字(它使用gmp作为依赖项。下面是一个示例：

library(Rmpfr)
1- (1/exp(mpfr(1000,precBits=1e5)))
[1] 0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999994924041102450543234708 ....