VGAM logit function overflows

Question

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?

Answer 1

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 ....