BUG: stable evaluation of expit

[argriffing]

closes #3385

[pv]

yield tests don't work in TestCase subclasses (nose weirdness...).

[argriffing]

Fixed the yield assertions. I didn't know they weren't allowed in TestCase, and it's disconcerting to me that they were just silent. I've probably committed some code in scipy or elsewhere with tests that are not doing anything...

[argriffing]

Is this the kind of bug that could be backported?

[pv]

Yes. However, @rgommers didn't branch 0.14.0 yet.

[coveralls]

Coverage remained the same when pulling 5a4fd30 on argriffing:expit-evaluation into 48b9d09 on scipy:master.

[coveralls]

Coverage remained the same when pulling 5a4fd30 on argriffing:expit-evaluation into 48b9d09 on scipy:master.

[argriffing]

I hope this fix is easy to review. Both conditional branches are mathematically the same, but they are numerically stable in different regimes (0<x vs x<0). The expit function is also called the logistic function (which is different from logit). See this blog post explaining how to evaluate it in a numerically stable way, branching according to the sign of x.

[WarrenWeckesser]

Looks good to me. There should also be tests for other data types, since the same problem occurs (at different thresholds) for float32 and float128:

In [50]: expit(np.array([88, 89], dtype=np.float32))
Out[50]: array([  1.,  nan], dtype=float32)

In [51]: expit(np.array([11356, 11357], dtype=np.float128))
Out[51]: array([ 1.0,  nan], dtype=float128)

[argriffing]

btw numpy and scipy are littered with bugs for float > 64 right? I never use these precisions, and the only time I see them is for reported bugs, so I assume that they are not truly supported beyond numpy.ndarray being able to physically store their representations.

[argriffing]

also my platform doesn't even have np.float96. Are all platforms guaranteed to have np.float128?

[argriffing]

@WarrenWeckesser added

[rgommers]

@pv correct. I planned to last night but there were too many PRs to merge the past two days. And at least I'd like to fix part of gh-3330; it makes little sense to produce a beta if I get a couple of hundred test errors myself.

[pv]

Thanks, LGTM

[larsmans]

I'm still surprised that it's evaluated this way. In scikit-learn, we had this implementation (in Cython), until I figured out that reducing expit to tanh was faster and stable enough (up to x=1000, at least):

def expit(x):
    z = .5 * x
    np.tanh(z, z)
    z += 1
    z *= .5
    return z

Has this ever been considered for the scipy.special implementation?