Python scipy.special.gammaincc() Examples

def longestrunones8(binin):
    ''' The focus of the test is the longest run of ones within M-bit blocks. The purpose of this test is to determine whether the length of the longest run of ones within the tested sequence is consistent with the length of the longest run of ones that would be expected in a random sequence. Note that an irregularity in the expected length of the longest run of ones implies that there is also an irregularity in the expected length of the longest run of zeroes. Long runs of zeroes were not evaluated separately due to a concern about statistical independence among the tests.'''
    m = 8
    k = 3
    pik = [0.2148, 0.3672, 0.2305, 0.1875]
    blocks = [binin[xs*m:m+xs*m:] for xs in range(len(binin) / m)]
    n = len(blocks)
    # append the string 01 to guarantee the length of 1
    counts1 = [xs+'01' for xs in blocks]
    counts = [xs.replace('0', ' ').split()
              for xs in counts1]  # split into all parts
    counts2 = [list(map(len, xx)) for xx in counts]
    counts4 = [(4 if xx > 4 else xx) for xx in map(max, counts2)]
    freqs = [counts4.count(spi) for spi in [1, 2, 3, 4]]
    chisqr1 = [(freqs[xx]-n*pik[xx])**2/(n*pik[xx]) for xx in range(4)]
    chisqr = reduce(su, chisqr1)
    pval = spc.gammaincc(k / 2.0, chisqr / 2.0)
    return pval 

def longestrunones128(binin):  # not well tested yet
    if len(binin) > 128:
        m = 128
        k = 5
        n = len(binin)
        pik = [0.1174, 0.2430, 0.2493, 0.1752, 0.1027, 0.1124]
        blocks = [binin[xs * m:m + xs * m:] for xs in range(len(binin) / m)]
        n = len(blocks)
        counts = [xs.replace('0', ' ').split() for xs in blocks]
        counts2 = [list(map(len, xx)) for xx in counts]
        counts3 = [(1 if xx < 1 else xx) for xx in map(max, counts2)]
        counts4 = [(4 if xx > 4 else xx) for xx in counts3]
        chisqr1 = [(counts4[xx] - n * pik[xx]) ** 2 / (n * pik[xx])
                   for xx in range(len(counts4))]
        chisqr = reduce(su, chisqr1)
        pval = spc.gammaincc(k / 2.0, chisqr / 2.0)
    else:
        print('longestrunones128 failed, too few bits:', len(binin))
        pval = 0
    return pval 

def longestrunones10000(binin):  # not well tested yet
    ''' The focus of the test is the longest run of ones within M-bit blocks. The purpose of this test is to determine whether the length of the longest run of ones within the tested sequence is consistent with the length of the longest run of ones that would be expected in a random sequence. Note that an irregularity in the expected length of the longest run of ones implies that there is also an irregularity in the expected length of the longest run of zeroes. Long runs of zeroes were not evaluated separately due to a concern about statistical independence among the tests.'''
    if len(binin) > 128:
        m = 10000
        k = 6
        pik = [0.0882, 0.2092, 0.2483, 0.1933, 0.1208, 0.0675, 0.0727]
        blocks = [binin[xs * m:m + xs * m:]
                  for xs in range(floor(len(binin) / m))]
        n = len(blocks)
        counts = [xs.replace('0', ' ').split() for xs in blocks]
        counts2 = [list(map(len, xx)) for xx in counts]
        counts3 = [(10 if xx < 10 else xx) for xx in map(max, counts2)]
        counts4 = [(16 if xx > 16 else xx) for xx in counts3]
        freqs = [counts4.count(spi) for spi in [10, 11, 12, 13, 14, 15, 16]]
        chisqr1 = [(freqs[xx] - n * pik[xx]) ** 2 / (n * pik[xx])
                   for xx in range(len(freqs))]
        chisqr = reduce(su, chisqr1)
        pval = spc.gammaincc(k / 2.0, chisqr / 2.0)
    else:
        print('longestrunones10000 failed, too few bits:', len(binin))
        pval = 0
    return pval

# test 2.06 

def randomexcursionstest(binin):
    ''' The focus of this test is the number of cycles having exactly K visits in a cumulative sum random walk. The cumulative sum random walk is found if partial sums of the (0,1) sequence are adjusted to (-1, +1). A random excursion of a random walk consists of a sequence of n steps of unit length taken at random that begin at and return to the origin. The purpose of this test is to determine if the number of visits to a state within a random walk exceeds what one would expect for a random sequence.'''
    xvals = [-4, -3, -2, -1, 1, 2, 3, 4]
    ss = [int(el) for el in binin]
    sc = list(map(sumi, ss))
    cumsum = np.cumsum(sc)
    cumsum = np.append(cumsum, 0)
    cumsum = np.append(0, cumsum)
    posi = np.where(cumsum == 0)[0]
    cycles = ([cumsum[posi[x]:posi[x+1]+1] for x in range(len(posi)-1)])
    j = len(cycles)
    sct = []
    for ii in cycles:
        sct.append(([len(np.where(ii == xx)[0]) for xx in xvals]))
    sct = np.transpose(np.clip(sct, 0, 5))
    su = []
    for ii in range(6):
        su.append([(xx == ii).sum() for xx in sct])
    su = np.transpose(su)
    pikt = ([([pik(uu, xx) for uu in range(6)]) for xx in xvals])
    # chitab=1.0*((su-j*pikt)**2)/(j*pikt)
    chitab = np.sum(1.0*(np.array(su)-j*np.array(pikt))
                    ** 2/(j*np.array(pikt)), axis=1)
    pval = ([spc.gammaincc(2.5, cs/2.0) for cs in chitab])
    return pval 

def inc_gamma(s, x):
    r"""The (upper) incomplete gamma function.

    Given by: :math:`\Gamma(s,x) = \int_x^{\infty} t^{s-1}\,e^{-t}\,{\rm d}t`

    Parameters
    ----------
    s : :class:`float`
        exponent in the integral
    x : :class:`numpy.ndarray`
        input values
    """
    if np.isclose(s, 0):
        return sps.exp1(x)
    if np.isclose(s, np.around(s)) and s < -0.5:
        return x ** (s - 1) * sps.expn(int(1 - np.around(s)), x)
    if s < 0:
        return (inc_gamma(s + 1, x) - x ** s * np.exp(-x)) / s
    return sps.gamma(s) * sps.gammaincc(s, x) 

def gammaincc(*args, **kwargs):
        from scipy.special import gammaincc
        return gammaincc(*args, **kwargs) 

def nonoverlappingtemplatematchingtest(binin, mat="000000001", num=9):
    ''' The focus of this test is the number of occurrences of pre-defined target substrings. The purpose of this test is to reject sequences that exhibit too many occurrences of a given non-periodic (aperiodic) pattern. For this test and for the Overlapping Template Matching test, an m-bit window is used to search for a specific m-bit pattern. If the pattern is not found, the window slides one bit position. For this test, when the pattern is found, the window is reset to the bit after the found pattern, and the search resumes.'''
    n = len(binin)
    m = len(mat)
    M = floor(n/num)
    blocks = [binin[xs*M:M+xs*M:] for xs in range(floor(n/M))]
    counts = [xx.count(mat) for xx in blocks]
    avg = 1.0 * (M-m+1)/2 ** m
    var = M*(2**-m - (2*m-1)*2**(-2*m))
    chisqr = reduce(su, [(xs - avg) ** 2 for xs in counts]) / var
    pval = spc.gammaincc(1.0 * len(blocks) / 2, chisqr / 2)
    return pval 

def serialtest(binin):
    m = int(log(len(binin), 2) - 3)
    ''' The focus of this test is the frequency of each and every overlapping m-bit pattern across the entire sequence. The purpose of this test is to determine whether the number of occurrences of the 2m m-bit overlapping patterns is approximately the same as would be expected for a random sequence. The pattern can overlap.'''
    n = len(binin)
    hbin = binin+binin[0:m-1:]
    f1a = [hbin[xs:m+xs:] for xs in range(n)]
    oo = set(f1a)
    f1 = [f1a.count(xs)**2 for xs in oo]
    f1 = list(map(f1a.count, oo))
    cou = f1a.count
    f2a = [hbin[xs:m-1+xs:] for xs in range(n)]
    f2 = [f2a.count(xs)**2 for xs in set(f2a)]
    f3a = [hbin[xs:m-2+xs:] for xs in range(n)]
    f3 = [f3a.count(xs)**2 for xs in set(f3a)]
    psim1 = 0
    psim2 = 0
    psim3 = 0
    if m >= 0:
        suss = reduce(su, f1)
        psim1 = 1.0 * 2 ** m * suss / n - n
    if m >= 1:
        suss = reduce(su, f2)
        psim2 = 1.0 * 2 ** (m - 1) * suss / n - n
    if m >= 2:
        suss = reduce(su, f3)
        psim3 = 1.0 * 2 ** (m - 2) * suss / n - n
    d1 = psim1-psim2
    d2 = psim1-2 * psim2 + psim3
    pval1 = spc.gammaincc(2 ** (m - 2), d1 / 2.0)
    pval2 = spc.gammaincc(2 ** (m - 3), d2 / 2.0)
    return [pval1, pval2] 

def aproximateentropytest(binin, m=5):
    ''' The focus of this test is the frequency of each and every overlapping m-bit pattern. The purpose of the test is to compare the frequency of overlapping blocks of two consecutive/adjacent lengths (m and m+1) against the expected result for a random sequence.'''
    n = len(binin)
    f1a = [(binin + binin[0:m - 1:])[xs:m + xs:] for xs in range(n)]
    f1 = [[xs, f1a.count(xs)] for xs in sorted(set(f1a))]
    f2a = [(binin + binin[0:m:])[xs:m + 1 + xs:] for xs in range(n)]
    f2 = [[xs, f2a.count(xs)] for xs in sorted(set(f2a))]
    c1 = [1.0 * f1[xs][1] / n for xs in range(len(f1))]
    c2 = [1.0 * f2[xs][1] / n for xs in range(len(f2))]
    phi1 = reduce(su, list(map(logo, c1)))
    phi2 = reduce(su, list(map(logo, c2)))
    apen = phi1 - phi2
    chisqr = 2.0 * n * (np.log(2) - apen)
    pval = spc.gammaincc(2 ** (m - 1), chisqr / 2.0)
    return pval 

def linearcomplexitytest(binin, m=500):
    ''' The focus of this test is the length of a generating feedback register. The purpose of this test is to determine whether or not the sequence is complex enough to be considered random. Random sequences are characterized by a longer feedback register. A short feedback register implies non-randomness.'''
    k = 6
    pi = [0.01047, 0.03125, 0.125, 0.5, 0.25, 0.0625, 0.020833]
    avg = 0.5*m + (1.0/36)*(9 + (-1)**(m + 1)) - (m/3.0 + 2.0/9)/2**m
    blocks = stringpart(binin, m)
    bign = len(blocks)
    lc = ([lincomplex(chunk) for chunk in blocks])
    t = ([-1.0*(((-1)**m)*(chunk-avg)+2.0/9) for chunk in lc])
    vg = np.histogram(t, bins=[-9999999999, -2.5, -
                               1.5, -0.5, 0.5, 1.5, 2.5, 9999999999])[0][::-1]
    im = ([((vg[ii]-bign*pi[ii])**2)/(bign*pi[ii]) for ii in range(7)])
    chisqr = reduce(su, im)
    pval = spc.gammaincc(k/2.0, chisqr/2.0)
    return pval 

def blockfrequencytest(binin, nu=20):
    ''' The focus of the test is the proportion of zeroes and ones within M-bit blocks. The purpose of this test is to determine whether the frequency of ones is an M-bit block is approximately M/2.'''
    ss = [int(el) for el in binin]
    tt = [1.0 * sum(ss[xs * nu:nu + xs * nu:]) /
          nu for xs in range(floor(len(ss) / nu))]
    uu = list(map(sus, tt))
    chisqr = 4 * nu * reduce(su, uu, 0)
    pval = spc.gammaincc(len(tt) / 2.0, chisqr / 2.0)
    return pval 

def erf(*args, **kwargs):
            from scipy.special import erf
            return erf(*args, **kwargs)


#    from scipy.special import lambertw, ellipe, gammaincc, gamma # fluids
#    from scipy.special import i1, i0, k1, k0, iv # ht
#    from scipy.special import hyp2f1    
#    if erf is None:
#        from scipy.special import erf 

def gammaincc(a, x):
        import mpmath
        return mpmath.gammainc(a, a=x, regularized=True) 

def gammaincc(a, x):
        import mpmath
        return mpmath.gammainc(a, a=x, regularized=True) 

def _sf(self, x, a):
        return sc.gammaincc(a, x) 

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

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 _cdf(self, x, a):
        return sc.gammaincc(a, 1.0 / x) 

def _cdf(self, x, beta):
        c = 0.5 * np.sign(x)
        # evaluating (.5 + c) first prevents numerical cancellation
        return (0.5 + c) - c * sc.gammaincc(1.0/beta, abs(x)**beta) 

def ggrb_int_cpl(a, Ec, Emin, Emax):

        # Gammaincc does not support quantities
        i1 = gammaincc(2 + a, old_div(Emin, Ec)) * gamma(2 + a)
        i2 = gammaincc(2 + a, old_div(Emax, Ec)) * gamma(2 + a)

        return -Ec * Ec * (i2 - i1) 

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

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

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 _cdf(self, x, a):
        return sc.gammaincc(a, 1.0 / x) 

def _cdf(self, x, beta):
        c = 0.5 * np.sign(x)
        # evaluating (.5 + c) first prevents numerical cancellation
        return (0.5 + c) - c * sc.gammaincc(1.0/beta, abs(x)**beta) 

def complete_gamma(a, z):
    """
    return 'complete' gamma function
    """
    return exp1(z) if a == 0 else gamma(a)*gammaincc(a, z) 

def test_gammaincc(self):
        assert_equal(cephes.gammaincc(5,0),1.0) 

def test_gammaincc(self):
        gicc = special.gammaincc(.5,.5)
        greal = 1 - special.gammainc(.5,.5)
        assert_almost_equal(gicc,greal,8) 

def test_gammainccnan(self):
        gama = special.gammaincc(-1,1)
        assert_(isnan(gama)) 

def _sf(self, x, a):
        return sc.gammaincc(a, x)