scipy/scipy — pull request #3545 by WarrenWeckesser (2014-04-13 19:09:28+00:00)
Title: BUG: stats: stats.levy.cdf with small arguments loses precision.

Change the implementation of stats.levy._cdf from

def _cdf(self, x):
return 2*(1-norm._cdf(1/sqrt(x)))

to

def _cdf(self, x):
return special.erfc(np.sqrt(0.5 / x))

The old implementation suffers from catastrophic cancellation when
x is small, because then norm._cdf(1/sqrt(x)) is approximately 1.

For example, with the old implementation:

In [56]: x = np.array([0.05, 0.015, 0.0125, 0.011, 0.01, 0.001])

In [57]: stats.levy.cdf(x)
Out[57]:
array([ 7.74421643e-06, 2.22044605e-16, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00])

The new implementation is:

In [58]: special.erfc(np.sqrt(0.5/x))
Out[58]:
array([ 7.74421643e-006, 3.21526273e-016, 3.74409738e-019,
1.50426951e-021, 1.52397060e-023, 1.79583278e-219])