Function gammai

pattern/metrics.py:830–879

Returns the incomplete gamma function for LOWER or UPPER tail.

(a, x, tail=UPPER)

    return exp(gammaln(x))

def gammai(a, x, tail=UPPER):
    """ Returns the incomplete gamma function for LOWER or UPPER tail.
    """
    
    # Series approximation to the incomplete gamma function.
    def gs(a, x, epsilon=3.e-7, iterations=700):
        ln = gammaln(a)
        s = 1.0 / a
        d = 1.0 / a
        for i in range(1, iterations):
            d = d * x / (a + i)
            s = s + d
            if abs(d) < abs(s) * epsilon:
                return (s * exp(-x + a * log(x) - ln), ln)
        raise StopIteration, (abs(d), abs(s) * epsilon)
    
    # Continued fraction approximation of the incomplete gamma function.
    def gf(a, x, epsilon=3.e-7, iterations=200):
        ln = gammaln(a)
        g0 = 0.0
        a0 = 1.0
        b0 = 0.0
        a1 = x
        b1 = 1.0
        f  = 1.0
        for i in range(1, iterations):
            a0 = (a1 + a0 * (i - a)) * f
            b0 = (b1 + b0 * (i - a)) * f
            a1 = x * a0 + a1 * i * f
            b1 = x * b0 + b1 * i * f
            if a1 != 0.0:
                f = 1.0 / a1
                g = b1 * f
                if abs((g - g0) / g) < epsilon:
                    return (g * exp(-x + a * log(x) - ln), ln)
                g0 = g
        raise StopIteration, (abs((g-g0) / g))

    if a <= 0.0:
        return 1.0
    if x <= 0.0:
        return 1.0
    if x < a + 1:
        if tail == LOWER:
            return gs(a, x)[0]
        return 1 - gs(a, x)[0]
    else:
        if tail == UPPER:
            return gf(a, x)[0]
        return 1 - gf(a, x)[0]

#--- COMPLEMENTARY ERROR FUNCTION ------------------------------------------------------------------
# Based on: http://www.johnkerl.org/python/sp_funcs_m.py.txt, Tom Loredo