logit and inverse logit functions for extreme values

Question

I need logit and inverse logit functions so that logit(inv_logit(n)) == n. I use numpy and here is what I have:

import numpy as np
def logit(p):
    return np.log(p) - np.log(1 - p)

def inv_logit(p):
    return np.exp(p) / (1 + np.exp(p))

And here are the values:

print logit(inv_logit(2)) 
2.0 

print logit(inv_logit(10))
10.0 

print logit(inv_logit(20))
20.000000018 #well, pretty close

print logit(inv_logit(50))
Warning: divide by zero encountered in log
inf 

Now let's test negative numbers

print logit(inv_logit(-10))
-10.0 
print logit(inv_logit(-20))
-20.0 
print logit(inv_logit(-200))
-200.0 
print logit(inv_logit(-500))
-500.0 
print logit(inv_logit(-2000))
Warning: divide by zero encountered in log
-inf 

So my questions are: what is the proper way to implement these functions so that the requirement logit(inv_logit(n)) == n will hold for any n in as wide a range as possible (at least [-1e4; 1e4)?

And also (and I'm sure this is connected to the first one), why are my function more stable with negative values, compared to the positive ones?

Answer 1

Nowadays, scipy has logit and expit (inverse logit) functions, eg

>>> from scipy.special import logit, expit
>>> import numpy as np
>>> logit([0, 0.25, 0.5, 0.75, 1])
array([       -inf, -1.09861229,  0.        ,  1.09861229,         inf])
>>> expit([-np.inf, -1.5, 0, 1.5, np.inf])
array([ 0.        ,  0.18242552,  0.5       ,  0.81757448,  1.        ])

Answer 2

You're running up against the precision limits for a IEEE 754 double-precision float. You'll need to use higher-precision numbers and operations if you want a larger range and a more precise domain.

>>> 1 + np.exp(-37)
1.0
>>> 1 + decimal.Decimal(-37).exp()
Decimal('1.000000000000000085330476257')