Python scipy.special.gammainc() Examples

def test_digamma_boundary():
    # Check that there isn't a jump in accuracy when we switch from
    # using the asymptotic series to the reflection formula.

    x = -np.logspace(300, -30, 100)
    y = np.array([-6.1, -5.9, 5.9, 6.1])
    x, y = np.meshgrid(x, y)
    z = (x + 1j*y).flatten()

    dataset = []
    with mpmath.workdps(30):
        for z0 in z:
            res = mpmath.digamma(z0)
            dataset.append((z0, complex(res)))
    dataset = np.asarray(dataset)

    FuncData(sc.digamma, dataset, 0, 1, rtol=1e-13).check()


# ------------------------------------------------------------------------------
# gammainc
# ------------------------------------------------------------------------------ 

def _upper_incomplete_gamma(z, a):
    """
    An implementation of the non-regularised upper incomplete gamma
    function. Computed using the relationship with the regularised 
    lower incomplete gamma function (scipy.special.gammainc). 
    Uses the recurrence relation wherever z<0.
    """
    n = int(-np.floor(z))
    if n > 0:
        z = z + n
        u_gamma = (1. - gammainc(z, a))*gamma(z)
        
        for i in range(n):
            z = z - 1.
            u_gamma = (u_gamma - np.power(a, z)*np.exp(-a))/z
        return u_gamma
    else:
        return (1. - gammainc(z, a))*gamma(z) 

def test_line():
    # Test on the line a = x where a simpler asymptotic expansion
    # (analog of DLMF 8.12.15) is available.

    def gammainc_line(x):
        c = np.array([-1/3, -1/540, 25/6048, 101/155520,
                      -3184811/3695155200, -2745493/8151736420])
        res = 0
        xfac = 1
        for ck in c:
            res -= ck*xfac
            xfac /= x
        res /= np.sqrt(2*np.pi*x)
        res += 0.5
        return res

    x = np.logspace(np.log10(25), 300, 500)
    a = x.copy()
    dataset = np.vstack((a, x, gammainc_line(x))).T

    FuncData(sc.gammainc, dataset, (0, 1), 2, rtol=1e-11).check() 

def mean(t, a, b):
        # TODO this is not tested yet.
        # tests:
        #    cemean(0., a, b)==mean(a, b, p)
        #    mean(t, a, 1., p)==mean(0., a, 1., p) == a
        # conditional excess mean
        # E[Y|y>t]
        # (conditional mean age at failure)
        # http://reliabilityanalyticstoolkit.appspot.com/conditional_distribution
        from scipy.special import gamma
        from scipy.special import gammainc
        # Regularized lower gamma
        print('not tested')

        v = 1. + 1. / b
        gv = gamma(v)
        L = np.power((t + .0) / a, b)
        cemean = a * gv * np.exp(L) * (1 - gammainc(v, t / a) / gv)

        return cemean 

def alpha_abs(self, x, y, n_sersic, r_eff, k_eff, center_x=0, center_y=0):
        """

        :param x:
        :param y:
        :param n_sersic:
        :param r_eff:
        :param k_eff:
        :param center_x:
        :param center_y:
        :return:
        """
        n = n_sersic
        x_red = self._x_reduced(x, y, n_sersic, r_eff, center_x, center_y)
        b = self.b_n(n_sersic)
        a_eff = self._alpha_eff(r_eff, n_sersic, k_eff)
        alpha = 2. * a_eff * x_red ** (-n) * (special.gammainc(2 * n, b * x_red))
        return alpha 

def test_gammainc(self):
        assert_equal(cephes.gammainc(5,0),0.0) 

def _cdf(self, x, df):
        return sc.gammainc(.5*df, .5*x**2) 

def test_gammainc(self):
        gama = special.gammainc(.5,.5)
        assert_almost_equal(gama,.7,1) 

def test_gammaincnan(self):
        gama = special.gammainc(-1,1)
        assert_(isnan(gama)) 

def test_gammainczero(self):
        # bad arg but zero integration limit
        gama = special.gammainc(-1,0)
        assert_equal(gama,0.0) 

def test_gammaincinf(self):
        gama = special.gammainc(0.5, np.inf)
        assert_equal(gama,1.0) 

def test_gammaincinv(self):
        y = special.gammaincinv(.4,.4)
        x = special.gammainc(.4,y)
        assert_almost_equal(x,0.4,1)
        y = special.gammainc(10, 0.05)
        x = special.gammaincinv(10, 2.5715803516000736e-20)
        assert_almost_equal(0.05, x, decimal=10)
        assert_almost_equal(y, 2.5715803516000736e-20, decimal=10)
        x = special.gammaincinv(50, 8.20754777388471303050299243573393e-18)
        assert_almost_equal(11.0, x, decimal=10) 

def test_975(self):
        # Regression test for ticket #975 -- switch point in algorithm
        # check that things work OK at the point, immediately next floats
        # around it, and a bit further away
        pts = [0.25,
               np.nextafter(0.25, 0), 0.25 - 1e-12,
               np.nextafter(0.25, 1), 0.25 + 1e-12]
        for xp in pts:
            y = special.gammaincinv(.4, xp)
            x = special.gammainc(0.4, y)
            assert_allclose(x, xp, rtol=1e-12) 

def test_gammainc_roundtrip():
    a = np.logspace(-5, 10, 100)
    x = np.logspace(-5, 10, 100)

    y = sc.gammaincinv(a, sc.gammainc(a, x))
    assert_allclose(x, y, rtol=1e-10) 

def test_gammainc(self):
        # Larger arguments are tested in test_data.py:test_local
        assert_mpmath_equal(sc.gammainc,
                            lambda z, b: mpmath.gammainc(z, b=b, regularized=True),
                            [Arg(0, 1e4, inclusive_a=False), Arg(0, 1e4)],
                            nan_ok=False, rtol=1e-11) 

def test_gammaincc(self):
        # Larger arguments are tested in test_data.py:test_local
        assert_mpmath_equal(sc.gammaincc,
                            lambda z, a: mpmath.gammainc(z, a=a, regularized=True),
                            [Arg(0, 1e4, inclusive_a=False), Arg(0, 1e4)],
                            nan_ok=False, rtol=1e-11) 

def _cdf(self, x, df):
        return sc.gammainc(.5*df, .5*x**2) 

def _cdf(self, x, a):
        fac = 0.5*sc.gammainc(a, abs(x))
        return np.where(x > 0, 0.5 + fac, 0.5 - fac) 

def _sf(self, x, a):
        fac = 0.5*sc.gammainc(a, abs(x))
        return np.where(x > 0, 0.5-fac, 0.5+fac) 

def _cdf(self, x, a):
        return sc.gammainc(a, x) 

def _sf(self, x, a, c):
        xc = x**c
        val1 = sc.gammainc(a, xc)
        val2 = sc.gammaincc(a, xc)
        return np.where(c > 0, val2, val1) 

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

def _cdf(self, x, c):
        return sc.gammainc(c, np.exp(x)) 

def _cdf(self, x):
        return sc.gammainc(1.5, x*x/2.0) 

def _cdf(self, x, nu):
        return sc.gammainc(nu, nu*x*x) 

def _cdf(self, x, a):
        return special.gammainc(a, x) 

def _cdf(self, x, c):
        return special.gammainc(c, numpy.exp(x)) 

def _cdf(self, x, a):
        fac = 0.5*special.gammainc(a,abs(x))
        return numpy.where(x>0,0.5+fac,0.5-fac) 

def _cdf(self, x, nu):
        return special.gammainc(nu,nu*x*x) 

def _cdf(self, x, a, c):
        val = special.gammainc(a, x**c)
        cond = c + 0*val
        return numpy.where(cond > 0, val, 1-val)