Comment gérer correctement Infs dans une fonction statistique?

8

Supposons que j'ai une fonction telle que:

f <- function(x){
  exp(x) / (1 + exp(x))
}

il est censé fonctionner pour toute valeur réelle de x, mais en fait, il renvoie NaN lorsque x est 710 ou plus. Je me demande quelle est la bonne façon de gérer ce problème. Je me rends compte qu'il est facile de le faire simplement renvoyer 1, mais ce n'est peut-être pas un bon comportement du point de vue d'un statisticien. Quelqu'un a-t-il des commentaires ou des suggestions?

r function numerics David Z
la source
Je ne sais pas si je pourrais faire confiance à des estimations de paramètres basées sur un modèle avec des valeurs d'influence aussi élevées dans la fonction. Vous pouvez vous attendre à ce que vos algorithmes de Newton-Raphson standard vous donnent des estimations de paramètres absurdes avec des valeurs de comme prédicteur linéaire dans les modèles de régression logistique. Les rapports de cotes peuvent être signalés comme une valeur infinie. De plus, je pense que vous pouvez inverser le test de score pour obtenir un intervalle de confiance valide pour le rapport de cotes. x
AdamO
Cela dépend vraiment du but vers lequel les valeurs sont tournées. pour les grands passe à ; cela peut être utile à certaines fins et pas très bon pour d'autres. exp(x)/(1+exp(x))x1−exp(−x)
Glen_b -Reinstate Monica

Réponses:

11

Dans ce cas, le NaN(pas un nombre) est renvoyé car le calcul des débordements exponentiels en arithmétique double précision.

Une expression algébriquement équivalente, développée dans une série MacLaurin autour de , est0

exp⁡(x)1+exp⁡(x)=11+exp⁡(−x)=1−exp⁡(−x)+exp⁡(−2x)−⋯.

Comme il s'agit d'une série alternée, l'erreur commise lors de la suppression d'un terme n'est pas supérieure à la taille du terme suivant. Ainsi, lorsque , l'erreur n'est pas supérieure à rapport à la valeur réelle. C'est beaucoup plus précis que n'importe quel calcul statistique doit être, donc vous êtes bien en remplaçant la valeur de retour par dans cette situation.x>710exp⁡(−710)≈10−308≈2−1024 1

Fait intéressant, Rne produira pas de NaNmoment où l'exponentielle se déverse . Ainsi, vous pouvez simplement choisir la version la plus fiable du calcul, en fonction du signe de x, comme dans

f <- function(x) ifelse(x < 0, exp(x) / (1 + exp(x)), 1 / (1 + exp(-x)))

Ce problème apparaît dans presque toutes les plates-formes informatiques (je n'ai pas encore vu d'exception) et elles varient dans la façon dont elles gèrent les débordements et les débordements. Les exponentielles sont connues pour créer ce genre de problèmes, mais elles ne sont pas seules. Par conséquent, il ne suffit pas d'avoir une solution R: un bon statisticien comprend les principes de l'arithmétique informatique et sait comment les utiliser pour détecter et contourner les particularités de son environnement informatique.

whuber
la source
1
Il peut être utile de souligner que lorsque environ, sera évalué à ( exactement ) en raison de l'arrondi en virgule flottante. De même, lorsque , évalué à , d'où le quotient produit une valeur exacte de . Les problèmes de précision lorsque sont astronomiquement plus petits! x<−361+exp⁡(x)1x>361+exp⁡(x)exp⁡(x)1|x|>710
whuber
1

D'autres ont déjà discuté des problèmes de calcul, je vais donc leur laisser cela. Puisque je suppose que vous travaillez avec R, je pensais que je soulignerais que le paquet de démarrage est livré avec sa propre fonction logit inverse à utiliser, qui est assez stable sur le plan des calculs:

require(boot) inv.logit(710)

semble évaluer à 1 comme souhaité.

Samuel Benidt
la source
1
Ou si vous souhaitez éviter d'introduire une dépendance de package, vous obtenez plogis(710)le même résultat. (En effet, inv.logitc'est juste un alias pour plogis.)
orizon