Переполнение функции логита VGAM

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

Обновление: это ошибка, несмотря на проблему с плавающей запятой, она должна возвращать 1. Действительно, plogis В этом случае следует использовать функцию, так как она правильно решает проблему. Автор и сопровождающий указаны в файле DESCRIPTION как Thomas Yee, вы должны отправить ему электронное письмо.

---

Ваша машина не может представлять плавающие точки, которые являются такими маленькими. Рассмотрим обратную функцию 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 ....