False CDF values for skew normal distribution

[pphilippos]

When x is sufficiently large, skewnorm.cdf outputs 0 instead of 1.

Reproducing code example:

A short code example that reproduces the problem:

from scipy.stats import skewnorm as sk2
import numpy as np
for x in np.linspace(-10,50,1000):
	print str(x) + " " + str(sk2.cdf(x,-1))

Error message:

No error messages. See the following graph for the issue:

Scipy/Numpy/Python version information:

Python 2.7.13
Scipy 0.19.1
Numpy 1.12.1

[pv]

Correct. The technical problem is that numerical integration of the pdf starts to fail once the support becomes too narrow compared to the interval. Can probably be fixed by just returning 1.0 when deep in the tail of the cdf.

[ev-br]

ISTM the best fix would be to finish off gh-7120, which implements the Owen's T function, and add the explicit form of the _cdf. This seems to fix the loss of precision, cf https://github.com/ev-br/scipy/tree/pr/7120. The only relevant commit on top of gh-7120 is 86b43ba

[ev-br]

Per the discussion in gh-8473: the naive formula contains the difference 2\Phi(x) - T(x,a), which suffers from the loss of precision for some parameters. A proper fix likely involves coming up with the strategy of computing the full expression, instead of subtracting two terms. Meanwhile, gh-8501 has a workaround.