Python math.gamma() Examples

def mvt_pdf(X, phi, K, nu):
    """
    Multivariate student-t density:
    output:
        the density of the given element
    input:
        x = parameter (d dimensional numpy array or scalar)
        mu = mean (d dimensional numpy array or scalar)
        K = scale matrix (dxd numpy array)
        nu = degrees of freedom
    """
    d = X.shape[-1]
    num = math.gamma((d + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * ((X - phi).dot(np.linalg.inv(K)).dot(np.transpose(X - phi))), -(d + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi, d / 2.) * pow(np.linalg.det(K), 0.5)
    return num / denom 

def expected_rs(n):
  """
  Calculates the expected (R/S)_n for white noise for a given n.

  This is used as a correction factor in the function hurst_rs. It uses the
  formula of Anis-Lloyd-Peters (see [h_3]_).

  Args:
    n (int):
      the value of n for which the expected (R/S)_n should be calculated

  Returns:
    float:
      expected (R/S)_n for white noise
  """
  front = (n - 0.5) / n
  i = np.arange(1,n)
  back = np.sum(np.sqrt((n - i) / i))
  if n <= 340:
    middle = math.gamma((n-1) * 0.5) / math.sqrt(math.pi) / math.gamma(n * 0.5)
  else:
    middle = 1.0 / math.sqrt(n * math.pi * 0.5)
  return front * middle * back 

def _levy_flight__(self, epoch, solution, prey):
        beta = 1
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        sigma_muy = np.power(gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)),1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy)
        v = np.random.normal(0, sigma_v)
        s = muy / np.power(np.abs(v), 1 / beta)
        # D is a random solution
        D = self._create_solution__(minmax=self.ID_MIN_PROBLEM)
        LB =0.001 * s * (solution[self.ID_POS] - prey[self.ID_POS])

        levy = D[self.ID_POS] * LB
        return levy

        # x_new = solution[0] + 1.0/np.sqrt(epoch+1) * np.sign(np.random.uniform() - 0.5) * levy
        # return x_new 

def AR_1_Coeff(data):
        """
        The autocovariance coefficient called as rho, for an AR(1) model can be calculated as shown here:

        .. math::

            \\rho(1) = \\frac{\\gamma(1)}{\\gamma(0)}

        For further information look for example in "Zeitreihenanalyse", pages 17, by Matti Schneider, Sebastian Mentemeier,
        SS 2010.

        :param data: numerical list
        :type data: list
        :return: autocorrelation coefficient
        :rtype: float
        """
        return TimeSeries.acf(data, 1) / TimeSeries.acf(data, 0) 

def _levy_flight__(self, solution, A, current_iter):
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        beta = 1
        sigma_muy = np.power(
            gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)),
            1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy)
        v = np.random.normal(0, sigma_v)
        s = muy / np.power(np.abs(v), 1 / beta)
        D = self._create_solution__(minmax=self.ID_MIN_PROBLEM)[self.ID_POS]
        LB = 0.01 * s * (solution - A)
        levy = D * LB
        # X_new = solution + 0.01*levy
        # X_new = solution + 1.0/np.sqrt(current_iter+1)*np.sign(np.random.random()-0.5)*levy
        return levy 

def _levy_flight__(self, epoch, solution, prey):
        beta = 1
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        sigma_muy = np.power(
            gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)),
            1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy**2)
        v = np.random.normal(0, sigma_v**2)
        s = muy / np.power(np.abs(v), 1 / beta)
        # D is a random solution
        D = self._create_solution__(minmax=self.ID_MAX_PROBLEM)
        LB = 0.01 * s * (solution[self.ID_POS] - prey[self.ID_POS])

        levy = D[self.ID_POS] * LB
        return levy

        # x_new = solution[0] + 1.0/np.sqrt(epoch+1) * np.sign(np.random.uniform() - 0.5) * levy
        # return x_new 

def _levy_flight__(self, epoch, solution, prey):
        beta = 1
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        sigma_muy = np.power(gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)), 1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy**2)
        v = np.random.normal(0, sigma_v**2)
        s = muy / np.power(np.abs(v), 1 / beta)
        # D is a random solution
        D = self._create_solution__(minmax=self.ID_MAX_PROBLEM)
        LB = 0.01 * s * (solution[self.ID_POS] - prey[self.ID_POS])

        levy = D[self.ID_POS] * LB
        return levy

        #x_new = solution[0] + 1.0/np.sqrt(epoch+1) * np.sign(np.random.uniform() - 0.5) * levy
        #return x_new 

def _levy_flight__(self, epoch, solution, prey):
        beta = 1
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        sigma_muy = np.power(gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)),1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy)
        v = np.random.normal(0, sigma_v)
        s = muy / np.power(np.abs(v), 1 / beta)
        # D is a random solution
        D = self._create_solution__()
        LB = 0.001 * s * (solution[self.ID_POS] - prey[self.ID_POS])

        levy = D[self.ID_POS] * LB
        return levy
        #x_new = solution[self.ID_POS] + 1.0/np.sqrt(epoch+1) * np.sign(np.random.uniform() - 0.5) * levy
        #return x_new 

def integrate(self, f, dot=numpy.dot):
        flt = numpy.vectorize(float)
        ref_vol = enr_volume(self.dim)
        return ref_vol * dot(f(flt(self.points).T), flt(self.weights))


# The closed formula is
#
#   2
#   * math.factorial(sum(alpha) + n - 1)
#   * prod([math.gamma((k + 1) / 2.0) for k in alpha])
#   / math.gamma((sum(alpha) + n) / 2).
#
# Care must be taken when evaluating this expression as numerator or denominator will
# quickly overflow. A better representation is via a recurrence. This is numerically
# stable and can easily be used for symbolic computation. 

def _levy_flight__(self, epoch, solution, prey):
        beta = 1
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        sigma_muy = np.power(gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)), 1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy)
        v = np.random.normal(0, sigma_v)
        s = muy / np.power(np.abs(v), 1 / beta)
        # D is a random solution
        D = self._create_solution__(minmax=self.ID_MAX_PROBLEM)
        LB = 0.001 * s * (solution[self.ID_POS] - prey[self.ID_POS])

        levy = D[self.ID_POS] * LB
        #return levy

        x_new = solution[0] + 1.0/np.sqrt(epoch+1) * np.sign(np.random.uniform() - 0.5) * levy
        return x_new 

def integrate_monomial_over_enr(k, symbolic=False):
    n = len(k)
    if any(a % 2 == 1 for a in k):
        return 0
    if all(a == 0 for a in k):
        return enr_volume(n, symbolic)

    # find first nonzero
    idx = next(i for i, j in enumerate(k) if j > 0)
    alpha = (k[idx] - 1) * (sum(k) + n - 1)
    k2 = k.copy()
    k2[idx] -= 2
    return integrate_monomial_over_enr(k2, symbolic) * alpha


# 2 * sqrt(pi) ** n * gamma(n) / gamma(frac(n, 2)) 

def _levy_flight__(self, epoch, solution, prey):
        beta = 1
        # muy and v are two random variables which follow normal distribution
        # sigma_muy : standard deviation of muy
        sigma_muy = np.power(gamma(1 + beta) * np.sin(np.pi * beta / 2) / (gamma((1 + beta) / 2) * beta * np.power(2, (beta - 1) / 2)), 1 / beta)
        # sigma_v : standard deviation of v
        sigma_v = 1
        muy = np.random.normal(0, sigma_muy)
        v = np.random.normal(0, sigma_v)
        s = muy / np.power(np.abs(v), 1 / beta)
        # D is a random solution
        D = self._create_solution__(minmax=self.ID_MAX_PROBLEM)
        LB = 0.01 * s * (solution[self.ID_POS] - prey[self.ID_POS])

        levy = D[self.ID_POS] * LB
        return levy

        #x_new = solution[0] + 1.0/np.sqrt(epoch+1) * np.sign(np.random.uniform() - 0.5) * levy
        #return x_new 

def crPlusMultifactor(N,M,wMat,p,c,aVec,alpha,rId):
    K = len(aVec)
    S = np.zeros([M,K])
    for k in range(0,K):        
        S[:,k] = np.random.gamma(aVec[k], 1/aVec[k], [M]) 
    W = wMat[rId,:]
    # Could replace tile with np.kron(W[:,0],np.ones([1,M])), but it's slow
    wS =  np.tile(W[:,0],[M,1]) + np.dot(S,np.transpose(W[:,1:]))
    pS = np.tile(p,[M,1])*wS
    H = np.random.poisson(pS,[M,N])
    lossIndicator = 1*np.greater_equal(H,1)
    lossDistribution = np.sort(np.dot(lossIndicator,c),axis=None)
    el,ul,var,es=util.computeRiskMeasures(M,lossDistribution,alpha)
    return el,ul,var,es          
      
# -----------------------
# Miscellaneous functions
# ----------------------- 

def test_integrate():
    moments = quadpy.tools.integrate(lambda x: [x ** k for k in range(5)], -1, +1)
    assert (moments == [2, 0, sympy.S(2) / 3, 0, sympy.S(2) / 5]).all()

    moments = quadpy.tools.integrate(
        lambda x: orthopy.line_segment.tree_legendre(x, 4, "monic", symbolic=True),
        -1,
        +1,
    )
    assert (moments == [2, 0, 0, 0, 0]).all()

    # Example from Gautschi's "How to and how not to" article
    moments = quadpy.tools.integrate(
        lambda x: [x ** k * sympy.exp(-(x ** 3) / 3) for k in range(5)], 0, sympy.oo
    )
    S = numpy.vectorize(sympy.S)
    gamma = numpy.vectorize(sympy.gamma)
    n = numpy.arange(5)
    reference = 3 ** (S(n - 2) / 3) * gamma(S(n + 1) / 3)
    assert numpy.all([sympy.simplify(m - r) == 0 for m, r in zip(moments, reference)]) 

def _sigma_r2_interp(self, r, a, gamma, rho0_r0_gamma, r_ani):
        """

        :param r:
        :param a:
        :param gamma:
        :param rho0_r0_gamma:
        :param r_ani:
        :return:
        """
        if not hasattr(self, '_interp_sigma_r2'):
            min_log = np.log10(self._min_integrate)
            max_log = np.log10(self._max_integrate)
            r_array = np.logspace(min_log, max_log, self._interp_grid_num)
            I_R_sigma2_array = []
            for r_i in r_array:
                I_R_sigma2_array.append(self._sigma_r2(r_i, a, gamma, rho0_r0_gamma, r_ani))
            self._interp_sigma_r2 = interp1d(np.log(r_array), np.array(I_R_sigma2_array), fill_value="extrapolate")
        return self._interp_sigma_r2(np.log(r)) 

def _read_out_params(self, kwargs_mass, kwargs_light, kwargs_anisotropy):
        """
        reads the relevant parameters out of the keyword arguments and transforms them to the conventions used in this
        class

        :param kwargs_mass: mass profile keyword arguments
        :param kwargs_light: light profile keyword arguments
        :param kwargs_anisotropy: anisotropy keyword arguments
        :return: a (Rs of Hernquist profile), gamma, rho0_r0_gamma, r_ani
        """
        if 'a' not in kwargs_light:
            kwargs_light['a'] = 0.551 * kwargs_light['r_eff']
        if 'rho0_r0_gamma' not in kwargs_mass:
            kwargs_mass['rho0_r0_gamma'] = self._rho0_r0_gamma(kwargs_mass['theta_E'], kwargs_mass['gamma'])
        a = kwargs_light['a']
        gamma = kwargs_mass['gamma']
        rho0_r0_gamma = kwargs_mass['rho0_r0_gamma']
        r_ani = kwargs_anisotropy['r_ani']
        return a, gamma, rho0_r0_gamma, r_ani 

def levy_flight(start: Tuple[float, float], alpha: float, param_lambda: float) -> Tuple[float, float]:
    """
    Perform a levy flight step.

    Arguments:
        start {Tuple[float, float]} -- The cuckoo's start position
        alpha {float} -- The step size
        param_lambda {float} -- lambda parameter of the levy distribution

    Returns:
        Tuple[float, float] -- The new position
    """

    def get_step_length():
        dividend = gamma(1 + param_lambda) * np.sin(np.pi * param_lambda / 2)
        divisor = gamma((1 + param_lambda) / 2) * param_lambda * np.power(2, (param_lambda - 1) / 2)
        sigma1 = np.power(dividend / divisor, 1 / param_lambda)

        sigma2 = 1

        u_vec = np.random.normal(0, sigma1, size=2)
        v_vec = np.random.normal(0, sigma2, size=2)

        step_length = u_vec / np.power(np.fabs(v_vec), 1 / param_lambda)

        return step_length

    return start + alpha * get_step_length() 

def _levy_flight_2__(self, solution=None, g_best=None):
        alpha = 0.01
        xichma_v = 1
        xichma_u = ((gamma(1 + 1.5) * np.sin(np.pi * 1.5 / 2)) / (gamma((1 + 1.5) / 2) * 1.5 * 2 ** ((1.5 - 1) / 2))) ** (1.0 / 1.5)
        levy_b = (np.random.normal(0, xichma_u ** 2)) / (np.sqrt(np.random.normal(0, xichma_v ** 2)) ** (1.0 / 1.5))
        return solution[self.ID_POS] + alpha * levy_b * (solution[self.ID_POS] - g_best[self.ID_POS]) 

def B(alpha: float, beta: float) -> float:
    """A normalizing constant so that the total probability is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 

def _levy_flight_HHO__(self, dimension=None, solution=None, prey=None):
        beta = 1.5
        sigma_muy = np.power((gamma(1+beta) * np.sin(np.pi*beta/2) / gamma((1+beta)/2)*beta*2**((beta-1)/2)) , 1.0/beta)
        u = np.random.uniform(self.domain_range[0], self.domain_range[1], dimension)
        v = np.random.uniform(self.domain_range[0], self.domain_range[1], dimension)

        LF_D =  0.01 * (sigma_muy * u) / np.power( np.abs(v), 1.0 / beta)
        LF_D2 = np.random.uniform(dimension) * LF_D
        Y = prey[self.ID_POS] + np.random.uniform(-1, 1) * np.abs( prey[self.ID_POS] - solution[self.ID_POS] )
        Z = Y + LF_D2
        return Z 

def poissonGammaAnalytic(N,c,a,b,alpha):
    pmfPoisson = np.zeros(N+1)
    q = np.divide(b,b+1)
    den = math.gamma(a)
    for k in range(0,N+1):
        num = np.divide(math.gamma(a+k),scipy.misc.factorial(k))
        pmfPoisson[k] = np.divide(num,den)*np.power(q,a)*np.power(1-q,k)  
    cdfPoisson = np.cumsum(pmfPoisson)
    varAnalytic = c*np.interp(alpha,cdfPoisson,np.linspace(0,N,N+1))
    esAnalytic = util.analyticExpectedShortfall(N,alpha,pmfPoisson,c)
    return pmfPoisson,cdfPoisson,varAnalytic,esAnalytic 

def gammaDensity(z,a,b):
    constant = np.divide(b**a,math.gamma(a))
    t1 = np.exp(-b*z)
    t2 = np.power(z,a-1)
    pdf =  constant*t1*t2
    return pdf 

def chi2Density(z,nu):
    g1 = math.gamma(nu/2)
    constant = np.multiply(g1,2**(nu/2))
    term1 = np.power(z,(nu/2)-1)
    term2 = np.exp(-z/2)   
    f = np.reciprocal(constant)*term1*term2
    return f 

def tDensity(z,mu,sigma,nu):
    g1 = math.gamma((nu+1)/2)
    g2 = math.gamma(nu/2)
    K = np.divide(g1,g2*np.sqrt(nu*math.pi)*sigma)
    power = np.divide((z-mu)**2,nu*(sigma**2))
    f = K*np.power(1+power,-(nu+1)/2)
    return f 

def computeBeta(a,b):
    Ga = math.gamma(a)
    Gb = math.gamma(b)
    Gab = math.gamma(a+b)
    return (Ga*Gb)/Gab 

def generateGamma(a,b,N):
    G1 = np.random.gamma(a,1,N)
    G2 = np.random.gamma(b,1,N)
    Z = np.divide(G1,G1+G2)  
    return Z 

def student_pdf_1d(X, mu, var, nu):
    if nu > 50:
        return gauss_pdf_1d(X, mu, var)
    num = math.gamma((1. + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * (1. / var * (X - mu) ** 2), -(1. + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * var, 0.5)
    return num / denom 

def student_pdf_1d(X, phi, K, nu):
    """
    Univariate student-t density:
    output:
        the density of the given element
    input:
        x = parameter scalar
        mu = mean
        var = scale matrix
        nu = degrees of freedom
    """
    num = math.gamma((1. + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * (1. / K * (X - phi) ** 2), -(1. + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * K, 0.5)
    return num / denom 

def B(alpha, beta):
    """a normalizing constant so that the total probability is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 

def crPlusOneFactor(N,M,w,p,c,v,alpha):
    S = np.random.gamma(v, 1/v, [M]) 
    wS =  np.transpose(np.tile(1-w + w*S,[N,1]))
    pS = np.tile(p,[M,1])*wS
    H = np.random.poisson(pS,[M,N])
    lossIndicator = 1*np.greater_equal(H,1)
    lossDistribution = np.sort(np.dot(lossIndicator,c),axis=None)
    el,ul,var,es=util.computeRiskMeasures(M,lossDistribution,alpha)
    return el,ul,var,es