Supongamos que tengo una función como:
f <- function(x){
exp(x) / (1 + exp(x))
}
se supone que funciona para cualquier valor real de x, pero en realidad devuelve NaN cuando x es 710 o mayor. Me pregunto cuál es la forma correcta de manejar este problema. Me doy cuenta de que es fácil hacer que solo devuelva 1, pero tal vez no sea un buen comportamiento desde el punto de vista de un estadístico. ¿Alguien tiene algunos comentarios o sugerencias?
Respuestas:
En este caso,
NaN
se devuelve (no un número) porque el cálculo de los desbordamientos exponenciales en aritmética de doble precisión.Una expresión algebraicamente equivalente, expandida en una serie de MacLaurin alrededor0 0 , es
Debido a que esta es una serie alterna, el error cometido al descartar cualquier término no es mayor que el tamaño del siguiente término. Así cuandox > 710 , el error no es mayor que Exp( - 710 ) ≈10- 308≈2- 1024 relativo al valor verdadero. Eso es mucho más preciso de lo que debe ser cualquier cálculo estadístico, por lo que está bien reemplazando el valor de retorno por1 en esta situación.
Curiosamente,
R
no producirá unNaN
cuando el exponencial se desborde . Por lo tanto, podría elegir la versión más confiable del cálculo, según el signo dex
, como enEste problema aparece en casi todas las plataformas informáticas (aún no he visto una excepción) y variarán en la forma en que manejan los desbordamientos y los desbordamientos. Los exponenciales son conocidos por crear este tipo de problemas, pero no están solos. Por lo tanto, no es suficiente tener una solución
R
: un buen estadístico comprende los principios de la aritmética informática y sabe cómo usarlos para detectar y solucionar las idiosincrasias de su entorno informático.fuente
Otros ya han discutido los problemas computacionales, por lo que se los dejaré Dado que supongo que está trabajando con R, pensé en señalar que el paquete de arranque viene con su propia función de logit inverso para su uso que es bastante estable desde el punto de vista computacional:
require(boot) inv.logit(710)
parece evaluar a 1 según lo deseado.
fuente
plogis(710)
logra el mismo resultado. (De hecho,inv.logit
es solo un alias paraplogis
.)