ENH: improve accuracy for ppf-cdf and sf-isf roundtrips for invgamma

[pbrod]

The accuracy of the invgamma ppf and cdf functions is poor for small x as exemplified here:

>>> import scipy.special as special
>>> import numpy as np 
>>> import scipy.stats as ss
>>> x = np.logspace(-2.6, 0)
>>> (ss.invgamma.ppf(ss.invgamma.cdf(x, 1), 1)-x)/x

 array([ -1.00000000e+00,  -1.00000000e+00,  -1.00000000e+00,
    -1.00000000e+00,  -1.00000000e+00,  -1.00000000e+00,
    -1.00000000e+00,  -1.00000000e+00,  -1.00000000e+00,
    -1.00000000e+00,  -1.00000000e+00,  -1.00000000e+00,
    -1.00000000e+00,  -1.00000000e+00,  -1.00000000e+00,
    -1.00000000e+00,  -1.00000000e+00,  -1.00000000e+00,
    -1.00000000e+00,  -1.00000000e+00,   1.05631197e-03,
     1.38355980e-05,  -9.84801496e-07,  -2.30789746e-09,
    -5.61840872e-10,  -4.08866497e-10,  -4.35294620e-11,
     2.29580766e-12,   1.48992764e-12,  -4.59710681e-13,
    -1.01250349e-13,  -3.55418986e-14,  -3.54414137e-15,
    -8.42943150e-15,  -3.12277800e-15,  -1.38181496e-15,
    -1.63052478e-15,   3.60749873e-16,  -4.25680444e-16,
     1.88361647e-16,   1.66698332e-16,   1.47526497e-16,
     3.91678795e-16,   2.31088097e-16,   2.04510849e-16,
     0.00000000e+00,   0.00000000e+00,   2.83506274e-16,
     0.00000000e+00,  -4.44089210e-16])

By replacing the _cdf and _ppf methods with

def _cdf(self, x):
    return special.gammaincc(a, 1.0/x)

def _ppf(self, q):
    return 1.0/special.gammainccinv(a, q)

one will obtain full machine accuracy as shown here:

>>> (ss.invgamma.ppf(ss.invgamma.cdf(x,1),1)-x)/x
array([  3.45302927e-16,  -1.52794964e-16,  -1.35222143e-16,
     2.39340716e-16,  -4.23628683e-16,   3.74907504e-16,
     1.65894853e-16,  -1.46815426e-16,   0.00000000e+00,
     0.00000000e+00,   2.03525116e-16,  -3.60235727e-16,
     0.00000000e+00,  -4.23209680e-16,   2.49691126e-16,
     0.00000000e+00,   3.91120568e-16,  -3.46138120e-16,
     0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    -6.36979987e-16,  -3.75814302e-16,   4.98888323e-16,
     2.94341064e-16,   1.30244573e-16,  -2.30530507e-16,
     0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
     0.00000000e+00,  -1.25147530e-16,   0.00000000e+00,
     1.96033291e-16,   0.00000000e+00,   0.00000000e+00,
     0.00000000e+00,   1.20249958e-16,   2.12840222e-16,
     0.00000000e+00,   0.00000000e+00,  -1.47526497e-16,
     2.61119197e-16,   0.00000000e+00,   2.04510849e-16,
    -1.80990228e-16,  -1.60174694e-16,   1.41753137e-16,
     0.00000000e+00,   0.00000000e+00])

The same applies to the _sf and _isf methods as shown here:

>>> x=np.logspace(2,100)
>>> (ss.invgamma.isf(ss.invgamma.sf(x, 1), 1)-x)/x 
array([ -5.40012479e-15,  -1.27329258e-14,   1.58618204e-11,
     1.07746929e-09,  -8.27903654e-08,   2.21222090e-05,
     7.99917193e-04,  -9.92800745e-02,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf,              inf,
                inf,              inf])

By replacing the _isf and _cdf methods with

def _isf(self, q):
    return 1.0/special.gammaincinv(a, q)

def _sf(self, x):
    return special.gammainc(a, 1.0/x)

one will get 1-3 digits accuracy for large x.

>>> (ss.invgamma.isf(ss.invgamma.sf(x, 1), 1)-x)/x 
array([ -5.68434189e-16,  -3.63797881e-16,   0.00000000e+00,
     2.23517418e-15,   1.52587891e-15,  -2.19726562e-15,
     9.37500000e-16,  -6.10215120e-03,  -4.23539248e-01,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03,  -4.48159239e-03,
    -4.48159239e-03,  -4.48159239e-03])

[person142]

In this case the round-trip accuracy could be misleading since gammaincinv and gammainccinv invert gammainc and gammaincc using Newton's method, i.e. they will probably be consistent but could be consistently wrong. (In fact gammainc and gammaincc both have some known weaknesses but will hopefully have less of them in a couple of days.)

It's still probably better than getting inf.