Skewnorm CDF inaccurate in the left tail

[dschmitz89]

SciPy currently falls back to quadrature integration of the PDF in the left tail of the skewnorm CDF for positive $a$ as the boost function has issues for low values. There seem to be catastrophic cancellations resulting in CDF values of 0. See the python example below. This is probably not a super serious issue but still wanted to make you aware in case you would like to fix it on your end.

import numpy as np
from scipy import special
import matplotlib.pyplot as plt

a = 1
x = np.linspace(-50, 0, 100000)

plt.semilogy(x, special._ufuncs._skewnorm_cdf(x, 0.0, 1.0, a))
plt.show()

Which results in the following plot:

Zooming in shows

For reference what the parameters of _skewnorm_cdf mean here the wrapper code:

template<typename Real>
Real
skewnorm_cdf_wrap(const Real x, const Real l, const Real sc, const Real sh)
{
    if (std::isfinite(x)) {
        return boost::math::cdf(
            boost::math::skew_normal_distribution<Real, StatsPolicy>(l, sc, sh), x);
    }
    // -inf => 0, inf => 1
    return 1 - std::signbit(x);
}

[jzmaddock]

I'll need to re-read the papers to see if there's a better way, although interestingly for shape=1 the cdf appears to simply be the square of the normal cdf which would at least solve that specific case.

[jzmaddock]

Notes to self, since I'm about to disappear on holiday ;)

Method T6 from "FAST AND ACCURATE CALCULATION OF OWEN’S T-FUNCTION" is interesting, but fails to provide an error bound.

"Evaluating the CDF of the Skew Normal Distribution", Amsler,2020 provides some interesting approximations for the far tails.