VGAM version 0.93
> logit(1000, inverse=T)
[1] 0 # it should be 1
The problem is here:
exp(1000 - log1p(exp(1000)))
Here log1p(exp(1000))
becomes Inf
So the numerical method it uses doesn't handle large numbers, compared to plogis
in base which works correctly.
Is it worth a bug report and where can I submit it?
Update: This is a bug, despite the floating point problem, it should return 1. Indeed, the plogis
function should be used in this case, as it correctly handles the problem. The author and maintainer are reported in the DESCRIPTION file as Thomas Yee, you should send him an email.
Your machine cannot represent floating points that are that small. Consider the inverse logit function:
inv.logit<-function(x) exp(x)/(1+exp(x))
Even at 500, it is very very small:
inv.logit(-500)
# 7.124576e-218
On my own machine, this is getting close to the limit of what the machine can represent. You can find this value .Machine$double.xmin
.
.Machine$double.xmin
# [1] 2.225074e-308
If you are truly interested in the exact value here, you will have to transform the numbers to a scale that can be represented on your computer.
Indeed, for large numbers the problem does not change. To be clear of what you are asking the machine to represent (when you ask for the inverse logit of 1000), try using the gmp
package. You are asking for the complement of the reciprocal of this number:
library(gmp)
exp(1)^as.bigz(1000)
Big Integer ('bigz') :
[1] 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
The problem is taking the reciprocal of this number, which is going to be very small, and unrepresentable on your computer.
You can actually use the Rmpfr
package to calculate this number (it uses gmp
as a dependency. Here is an example:
library(Rmpfr)
1- (1/exp(mpfr(1000,precBits=1e5)))
[1] 0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999994924041102450543234708 ....