Python scipy.special.gamma() Examples

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 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 _pdf_skip(self, x, dfn, dfd, nc):
        # ncf.pdf(x, df1, df2, nc) = exp(nc/2 + nc*df1*x/(2*(df1*x+df2))) *
        #             df1**(df1/2) * df2**(df2/2) * x**(df1/2-1) *
        #             (df2+df1*x)**(-(df1+df2)/2) *
        #             gamma(df1/2)*gamma(1+df2/2) *
        #             L^{v1/2-1}^{v2/2}(-nc*v1*x/(2*(v1*x+v2))) /
        #             (B(v1/2, v2/2) * gamma((v1+v2)/2))
        n1, n2 = dfn, dfd
        term = -nc/2+nc*n1*x/(2*(n2+n1*x)) + sc.gammaln(n1/2.)+sc.gammaln(1+n2/2.)
        term -= sc.gammaln((n1+n2)/2.0)
        Px = np.exp(term)
        Px *= n1**(n1/2) * n2**(n2/2) * x**(n1/2-1)
        Px *= (n2+n1*x)**(-(n1+n2)/2)
        Px *= sc.assoc_laguerre(-nc*n1*x/(2.0*(n2+n1*x)), n2/2, n1/2-1)
        Px /= sc.beta(n1/2, n2/2)
        # This function does not have a return.  Drop it for now, the generic
        # function seems to work OK. 

def _pdf(self, x, df, nc):
        # nct.pdf(x, df, nc) =
        #                               df**(df/2) * gamma(df+1)
        #                ----------------------------------------------------
        #                2**df*exp(nc**2/2) * (df+x**2)**(df/2) * gamma(df/2)
        n = df*1.0
        nc = nc*1.0
        x2 = x*x
        ncx2 = nc*nc*x2
        fac1 = n + x2
        trm1 = n/2.*np.log(n) + sc.gammaln(n+1)
        trm1 -= n*np.log(2)+nc*nc/2.+(n/2.)*np.log(fac1)+sc.gammaln(n/2.)
        Px = np.exp(trm1)
        valF = ncx2 / (2*fac1)
        trm1 = np.sqrt(2)*nc*x*sc.hyp1f1(n/2+1, 1.5, valF)
        trm1 /= np.asarray(fac1*sc.gamma((n+1)/2))
        trm2 = sc.hyp1f1((n+1)/2, 0.5, valF)
        trm2 /= np.asarray(np.sqrt(fac1)*sc.gamma(n/2+1))
        Px *= trm1+trm2
        return Px 

def _munp(self, n, beta, m):
        """
        Returns the n-th non-central moment of the crystalball function.
        """
        N = 1.0 / (m/beta / (m-1) * np.exp(-beta**2 / 2.0) + _norm_pdf_C * _norm_cdf(beta))

        def n_th_moment(n, beta, m):
            """
            Returns n-th moment. Defined only if n+1 < m
            Function cannot broadcast due to the loop over n
            """
            A = (m/beta)**m * np.exp(-beta**2 / 2.0)
            B = m/beta - beta
            rhs = 2**((n-1)/2.0) * sc.gamma((n+1)/2) * (1.0 + (-1)**n * sc.gammainc((n+1)/2, beta**2 / 2))
            lhs = np.zeros(rhs.shape)
            for k in range(n + 1):
                lhs += sc.binom(n, k) * B**(n-k) * (-1)**k / (m - k - 1) * (m/beta)**(-m + k + 1)
            return A * lhs + rhs

        return N * _lazywhere(np.atleast_1d(n + 1 < m),
                              (n, beta, m),
                              np.vectorize(n_th_moment, otypes=[np.float]),
                              np.inf) 

def compute_finite_difference_coefficients(derivative_order, grid_size):

    # from http://www.scientificpython.net/pyblog/uniform-finite-differences-all-orders

    n = 2*grid_size -1
    A = np.tile(np.arange(grid_size), (n,1)).T
    B = np.tile(np.arange(1-grid_size,grid_size), (grid_size,1))
    M = (B**A)/gamma(A+1)

    r = np.zeros(grid_size)
    r[derivative_order] = 1

    D = np.zeros((grid_size, grid_size))
    for k in range(grid_size):
        indexes = k + np.arange(grid_size)
        D[:,k] = solve(M[:,indexes], r)

    return D

#################################################################################################### 

def cgmy_chf(u, t, c, g, m, y):
    return np.exp(c * t * gamma(-y) * ((m - 1j * u) ** y - m ** y + (g + 1j * u) ** y - g ** y)) 

def factorial(n):
    return gamma(n + 1) 

def volume_unit_ball(n_dims):
    return (pi ** (.5 * n_dims)) / gamma(.5 * n_dims + 1) 

def sigma_e(re, n, q, Ftot=1., c0=0.):
    """
    Surface brightess at the effective radius, re, given total flux
    """
    from scipy.special import gamma
    kap = get_kappa(n)
    return Ftot/(2*np.pi*re**2*np.exp(kap)*n*kap**(-2*n)*gamma(2*n)*q/Rc(c0)), kap 

def phase_covariance_von_karman(r0, L0):
	'''Return a Field generator for the phase covariance function for Von Karman turbulence.

	Parameters
	----------
	r0 : scalar
		The Fried parameter.
	L0 : scalar
		The outer scale.

	Returns
	-------
	Field generator
		The phase covariance Field generator.
	'''
	def func(grid):
		r = grid.as_('polar').r + 1e-10

		a = (L0 / r0)**(5 / 3)
		b = gamma(11 / 6) / (2**(5 / 6) * np.pi**(8 / 3))
		c = (24 / 5 * gamma(6 / 5))**(5 / 6)
		d = (2 * np.pi * r / L0)**(5 / 6)
		e = kv(5 / 6, 2 * np.pi * r / L0)

		return Field(a * b * c * d * e, grid)
	return func 

def apply(self, solution):
        corr = solution.corr
        err = solution.err
        m = solution.model
        cond = solution.conditions
        undec = solution.undec
        evolution = solution.evolution
        # Make gamma distribution
        gamma_mean = self.pausestop - self.pausestart
        gamma_var = pow(self.pauseblurwidth, 2)
        shape = gamma_mean**2/gamma_var
        scale = gamma_var/gamma_mean
        gamma_pdf = lambda t : 1/(sp_gamma(shape)*(scale**shape)) * t**(shape-1) * np.exp(-t/scale)
        gamma_vals = np.asarray([gamma_pdf(t) for t in m.t_domain() - self.pausestart if t >= 0])
        sumgamma = np.sum(gamma_vals)
        gamma_start = next(i for i,t in enumerate(m.t_domain() - self.pausestart) if t >= 0)

        # Generate first part of pdf (before the pause)
        newcorr = np.zeros(m.t_domain().shape, dtype=corr.dtype)
        newerr = np.zeros(m.t_domain().shape, dtype=err.dtype)
        # Generate pdf after the pause
        for i,t in enumerate(m.t_domain()):
            #print(np.sum(newcorr)+np.sum(newerr))
            if 0 <= t < self.pausestart:
                newcorr[i] = corr[i]
                newerr[i] = err[i]
            elif self.pausestart <= t:
                newcorr[i:] += corr[gamma_start:len(corr)-(i-gamma_start)]*gamma_vals[int(i-gamma_start)]/sumgamma
                newerr[i:] += err[gamma_start:len(corr)-(i-gamma_start)]*gamma_vals[int(i-gamma_start)]/sumgamma
            else:
                raise ValueError("Invalid domain")
        return Solution(newcorr, newerr, m, cond, undec, evolution) 

def __repr__(self):
        return "{0}(gamma={1}, metric={2})".format(
            self.__class__.__name__, self.gamma, self.metric) 

def hyperparameter_gamma(self):
        return Hyperparameter("gamma", "numeric", self.gamma_bounds) 

def phase_structure_function_von_karman(r0, L0):
	'''Return a Field generator for the phase structure function for Von Karman turbulence.

	Parameters
	----------
	r0 : scalar
		The Fried parameter.
	L0 : scalar
		The outer scale.

	Returns
	-------
	Field generator
		The phase structure Field generator.
	'''
	def func(grid):
		r = grid.as_('polar').r + 1e-10

		a = (L0 / r0)**(5 / 3)
		b = 2**(1 / 6) * gamma(11 / 6) / np.pi**(8 / 3)
		c = (24 / 5 * gamma(6 / 5))**(5 / 6)
		d = gamma(5 / 6) / 2**(1 / 6)
		e = (2 * np.pi * r / L0)**(5 / 6)
		f = kv(5 / 6, 2 * np.pi * r / L0)

		return Field(a * b * c * (d - e * f), grid)
	return func 

def test_gammaln(self):
        gamln = special.gammaln(3)
        lngam = log(special.gamma(3))
        assert_almost_equal(gamln,lngam,8) 

def test_beta(self):
        bet = special.beta(2,4)
        betg = (special.gamma(2)*special.gamma(4))/special.gamma(6)
        assert_almost_equal(bet,betg,8) 

def test_gamma(self):
        assert_equal(cephes.gamma(5),24.0) 

def fromspline(cls, xk, cvals, order, fill=0.0):
        N = len(xk)-1
        sivals = np.empty((order+1, N), dtype=float)
        for m in xrange(order, -1, -1):
            fact = spec.gamma(m+1)
            res = _fitpack._bspleval(xk[:-1], xk, cvals, order, m)
            res /= fact
            sivals[order-m, :] = res
        return cls(sivals, xk, fill=fill) 

def dogft(s: np.ndarray,
          w: np.ndarray,
          order: int,
          dt: int) -> np.ndarray:
    """
    Fourier transformed DOG function.

    Parameters
    ----------
    s : numpy.ndarray
        Scales.
    w : numpy.ndarray
        Angular frequencies.
    order : int
        Wavelet order.
    dt : int
        Time step.

    Returns
    -------
    numpy.ndarray
        Normalized Fourier transformed DOG function.
    """

    p = - (0.0 + 1.0j) ** order / np.sqrt(gamma(order + 0.5))
    wavelet = np.zeros((s.shape[0], w.shape[0]), dtype=np.complex128)

    for i in range(s.shape[0]):
        n = normalization(s[i], dt)
        h = s[i] * w
        wavelet[i] = n * p * h ** order * np.exp(-h ** 2 / 2.0)

    return wavelet 

def cgmy_mgf(u, t, c, g, m, y):
    return np.exp(c * t * gamma(-y) * ((m - u) ** y - m ** y + (g + u) ** y - g ** y)) 

def compute_features(img_norm):
    features = []
    alpha, beta_left, beta_right = estimate_aggd_params(img_norm)

    features.extend([ alpha, (beta_left+beta_right)/2 ])

    for x_shift, y_shift in ((0,1), (1,0), (1,1), (1,-1)):
        img_pair_products  = img_norm * numpy.roll(numpy.roll(img_norm, y_shift, axis=0), x_shift, axis=1)
        alpha, beta_left, beta_right = estimate_aggd_params(img_pair_products)
        eta = (beta_right - beta_left) * (gamma(2.0/alpha) / gamma(1.0/alpha))
        features.extend([ alpha, eta, beta_left, beta_right ])

    return features 

def mean(a, b):
    """Continuous mean. Theoretically at most 1 step below discretized mean

    `E[T ] <= E[Td] + 1` true for positive distributions.

    :param a: Alpha
    :param b: Beta
    :return: `a * gamma(1.0 + 1.0 / b)`
    """
    from scipy.special import gamma
    return a * gamma(1.0 + 1.0 / b) 

def mean(a, b):
    """ Continuous mean. at most 1 step below discretized mean 

    `E[T ] <= E[Td] + 1` true for positive distributions.
    """
    from scipy.special import gamma
    return a * gamma(1.0 + 1.0 / b) 

def _hs(k, cs, rho, omega):
    c0 = (cs * cs * (1 + rho * rho) / (1 - rho * rho) /
          (1 - 2 * rho * rho * cos(2 * omega) + rho ** 4))
    gamma = (1 - rho * rho) / (1 + rho * rho) / tan(omega)
    ak = abs(k)
    return c0 * rho ** ak * (cos(omega * ak) + gamma * sin(omega * ak)) 

def _chi2_moment(n, df):
    # (raw) moments of \chi^2(df)
    return (2**n) * gamma(n + df/2.) / gamma(df/2.) 

def bghfactor(df):
    return np.power(2.0, 1-df*0.5) / sps_gamma(df*0.5) 

def chi_pdf(x, df):
    tmp = (df-1.)*np_log(x) + (-x*x*0.5) - (df*0.5-1)*np_log(2.0) \
          - sps_gammaln(df*0.5)
    return np_exp(tmp)
    #return x**(df-1.)*np_exp(-x*x*0.5)/(2.0)**(df*0.5-1)/sps_gamma(df*0.5) 

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

def _dog_function(order: int,
                      x_in: float) -> float:
        """
        Returns
        -------
        float
            DOG function.
        """

        p_hpoly = hermite(order)[int(x_in / np.power(2, 0.5))]
        herm = p_hpoly / (np.power(2, float(order) / 2))

        return ((-1)**(order+1)) / np.sqrt(gamma(order + 0.5)) * herm