Python math.lgamma() Examples

def PDF(self, x):
        if (self.alpha == 0.0) or (self.beta == 0.0):
            if (x != 0.0):
                return 0.0
            else:
                return float("inf")
        elif (x < 0.0):
            return 0.0
        elif (self.alpha == 1.0) and (x == 0.0):
            return self.beta
        else:
            try:
                term1a = math.log(x/self.beta) * (self.alpha - 1.0)
            except ValueError:
                term1a = float("-inf") * (self.alpha - 1.0)
            term1 = term1a - math.log(self.beta)
            term2 = -x/self.beta
            term3 = math.lgamma(self.alpha)
            return math.exp(term1 + (term2 - term3)) 

def __init__(self, alpha, beta, name, errcodeBase, pos):
        self.name = name
        self.errcodeBase = errcodeBase
        self.pos = pos
        # first shape parameter
        self.alpha = alpha
        # second shape parameter
        self.beta  = beta
        if (self.alpha <= 0.0) or (self.beta <= 0.0):
            raise PFARuntimeException("invalid parameterization", self.errcodeBase + 0, self.name, self.pos)
        # normalization factor
        self.Z = math.lgamma(self.alpha) + math.lgamma(self.beta) \
                 - math.lgamma(self.alpha + self.beta)
        # tolerance
        self.epsilon = 1e-15
        # max Iterations
        self.maxIter = 1000 

def multinomLog2(selectors):
    """
    Function calculates logarithm 2 of a kind of multinom.

    selectors: list of integers
    """

    ln2 = 0.69314718055994528622
    noAll = sum(selectors)
    lgNf = math.lgamma(noAll + 1.0) / ln2  # log2(N!)

    lgnFac = []
    for selector in selectors:
        if selector == 0 or selector == 1:
            lgnFac.append(0.0)
        elif selector == 2:
            lgnFac.append(1.0)
        elif selector == noAll:
            lgnFac.append(lgNf)
        else:
            lgnFac.append(math.lgamma(selector + 1.0) / ln2)
    return lgNf - sum(lgnFac) 

def regularizedGammaP(a, x):
    if (x < 0.0) or (a <= 0.0):
        raise Exception
    maxIter = 1000
    epsilon = 1e-15
    n  = 0.0
    an = 1.0/a
    total = an
    while (abs(an/total) > epsilon) and \
          (n < maxIter) and \
          (total < float("inf")):
        n = n + 1.0
        an = an * (x / (a + n))
        total = total + an
    if n >= maxIter:
        raise Exception
    elif total == float("inf"):
        ans = 1.0
    else:
        ans = math.exp(-x + (a*math.log(x)) - math.lgamma(a)) * total
    return ans

# http://www.johnkerl.org/python/sp_funcs_m.py.txt 

def gser(a, x, itmax=700, eps=3.0e-7):
    """Series approximation to the incomplete gamma function."""
    gln = math.lgamma(a)
    if (x < 0.0):
        raise Exception
    if (x == 0.0):
        return [0.0]
    ap = a
    total = 1.0 / a
    delta = total
    n = 1
    while n <= itmax:
        ap = ap + 1.0
        delta = delta * x / ap
        total = total + delta
        if (abs(delta) < abs(total)*eps):
            return [total * math.exp(-x + a*math.log(x) - gln), gln]
        n = n + 1
    raise Exception

# http://www.johnkerl.org/python/sp_funcs_m.py.txt
# used by incompleteGamma 

def log_betabin(k, n, alpha, beta):
    try:
        c = math.lgamma(alpha + beta) - math.lgamma(alpha) - math.lgamma(beta)
        # logger.debug("c: %f" % c)
    except:
        print('alpha = {}, beta = {}'.format(alpha, beta))
        raise
    if isinstance(k, (list, np.ndarray)):
        if len(k) != len(n):
            print('length of k in %d and length of n is %d' % (len(k), len(n)))
            raise ValueError
        lbeta = []
        for ki, ni in zip(k, n):
            lbeta.append(math.lgamma(ki + alpha) + math.lgamma(ni - ki + beta) - math.lgamma(ni + alpha + beta) + c)
        return np.array(lbeta)
    else:
        return math.lgamma(k + alpha) + math.lgamma(max(1, n - k + beta)) - math.lgamma(n + alpha + beta) + c 

def nChooseK(n, k):
    # is n an integer?
    nInt = (math.floor(n) == n)
    if n == k or k == 0:
        return 1
    if (n < k) or (k < 0):
        raise Exception
    if (nInt) and (n < 0.0):
        b = pow(-1.0, k) * math.exp(math.lgamma(abs(n + k)) \
                                  - math.lgamma(k + 1.0)    \
                                  - math.lgamma(abs(n)))
        return round(b)
    if (n >= k):
        b = math.exp(math.lgamma(n + 1.0) - math.lgamma(k + 1.0) \
                   - math.lgamma(n - k + 1.0))
        return round(b)
    if not (nInt) and (n < k):
        b = (1.0/math.pi) * math.exp(math.lgamma(n + 1.0) \
                                   - math.lgamma(k + 1)   \
                                   + math.lgamma(k - n)   \
                   + math.log(math.sin(math.pi * (n - k + 1.0))))
        return round(b)
    return 0.0 

def log_betabin(k,n,alpha,beta):
    import math
    try:
        Const =  math.lgamma(alpha + beta) - math.lgamma(alpha) - math.lgamma(beta)
    except:
        print 'alpha = {}, beta = {}'.format(alpha,beta)
    if isinstance(k,list) or isinstance(k,np.ndarray):
        if len(k)!=len(n):
            print 'length of k is %d and length of n is %d'%(len(k),len(n))
            raise ValueError
        lbeta = []
        for ki,ni in zip(k,n):
            # lbeta.append(math.lgamma(ni+1)- math.lgamma(ki+1) - math.lgamma(ni-ki+1) + math.lgamma(ki+alpha) + math.lgamma(ni-ki+beta) - math.lgamma(ni+alpha+beta) + Const)
            lbeta.append(math.lgamma(ki+alpha) + math.lgamma(ni-ki+beta) - math.lgamma(ni+alpha+beta) + Const)
        return np.array(lbeta)
    else:
        return math.lgamma(k+alpha) + math.lgamma(n-k+beta) - math.lgamma(n+alpha+beta) + Const
        # return math.lgamma(n+1)- math.lgamma(k+1) - math.lgamma(n-k+1) + math.lgamma(k+alpha) + math.lgamma(n-k+beta) - math.lgamma(n+alpha+beta) + Const 

def _gser(a,x):
    "Series representation of Gamma. NumRec sect 6.1."
    ITMAX=100
    EPS=3.e-7

    gln=lgamma(a)
    assert(x>=0),'x < 0 in gser'
    if x == 0 : return 0,gln

    ap = a
    delt = sum = 1./a
    for i in range(ITMAX):
        ap=ap+1.
        delt=delt*x/ap
        sum=sum+delt
        if abs(delt) < abs(sum)*EPS: break
    else:
        print('a too large, ITMAX too small in gser')
    gamser=sum*np.exp(-x+a*np.log(x)-gln)
    return gamser,gln 

def testFloatEmpty(self):
    x = np.empty((2, 0, 5), dtype=np.float32)
    self._compareBoth(x, np.abs, tf.abs)
    self._compareBoth(x, np.abs, _ABS)
    self._compareBoth(x, np.negative, tf.neg)
    self._compareBoth(x, np.negative, _NEG)
    self._compareBoth(x, self._inv, tf.inv)
    self._compareBoth(x, np.square, tf.square)
    self._compareBoth(x, np.sqrt, tf.sqrt)
    self._compareBoth(x, self._rsqrt, tf.rsqrt)
    self._compareBoth(x, np.exp, tf.exp)
    self._compareBoth(x, np.log, tf.log)
    self._compareBoth(x, np.log1p, tf.log1p)
    self._compareBoth(x, np.tanh, tf.tanh)
    self._compareBoth(x, self._sigmoid, tf.sigmoid)
    self._compareBoth(x, np.sign, tf.sign)
    self._compareBoth(x, np.sin, tf.sin)
    self._compareBoth(x, np.cos, tf.cos)
    # Can't use vectorize below, so just use some arbitrary function
    self._compareBoth(x, np.sign, tf.lgamma)
    self._compareBoth(x, np.sign, tf.erf)
    self._compareBoth(x, np.sign, tf.erfc)
    self._compareBoth(x, np.tan, tf.tan)
    self._compareBoth(x, np.arcsin, tf.asin)
    self._compareBoth(x, np.arccos, tf.acos)
    self._compareBoth(x, np.arctan, tf.atan)

    self._compareBothSparse(x, np.abs, tf.abs)
    self._compareBothSparse(x, np.negative, tf.neg)
    self._compareBothSparse(x, np.square, tf.square)
    self._compareBothSparse(x, np.sqrt, tf.sqrt, tol=1e-3)
    self._compareBothSparse(x, np.tanh, tf.tanh)
    self._compareBothSparse(x, np.sign, tf.sign)
    self._compareBothSparse(x, np.sign, tf.erf) 

def test_gamma_lgamma(self):
        tolerance = 13
        for x in itertools.count(0.001, 0.001):
            if x > 5.0:
                break
            self.assertAlmostEqual(math.lgamma(x), math.log(math.gamma(x)), places=tolerance)
            self.assertAlmostEqual(math.lgamma(x*x), math.log(math.gamma(x*x)), places=tolerance)
            self.assertAlmostEqual(math.lgamma(2.0**x), math.log(math.gamma(2.0**x)), places=tolerance)

            # Test negative values too, but not integers
            if x % 1.0 != 0.0:
                self.assertAlmostEqual(math.lgamma(-x), math.log(abs(math.gamma(-x))), places=tolerance)
                self.assertAlmostEqual(math.lgamma(-x*x), math.log(abs(math.gamma(-x*x))), places=tolerance)
                self.assertAlmostEqual(math.lgamma(-2.0**x), math.log(abs(math.gamma(-2.0**x))), places=tolerance) 

def testDoubleBasic(self):
    x = np.arange(-3, 3).reshape(1, 3, 2).astype(np.float64)
    y = (x + .5).astype(np.float64)    # no zero
    z = (x + 15.5).astype(np.float64)  # all positive
    k = np.arange(-0.90, 0.90, 0.35).reshape(1, 3, 2).astype(np.float64) # between -1 and 1
    self._compareBoth(x, np.abs, tf.abs)
    self._compareBoth(x, np.abs, _ABS)
    self._compareBoth(x, np.negative, tf.neg)
    self._compareBoth(x, np.negative, _NEG)
    self._compareBoth(y, self._inv, tf.inv)
    self._compareBoth(x, np.square, tf.square)
    self._compareBoth(z, np.sqrt, tf.sqrt)
    self._compareBoth(z, self._rsqrt, tf.rsqrt)
    self._compareBoth(x, np.exp, tf.exp)
    self._compareBoth(z, np.log, tf.log)
    self._compareBoth(z, np.log1p, tf.log1p)
    self._compareBoth(x, np.tanh, tf.tanh)
    self._compareBoth(x, self._sigmoid, tf.sigmoid)
    self._compareBoth(y, np.sign, tf.sign)
    self._compareBoth(x, np.sin, tf.sin)
    self._compareBoth(x, np.cos, tf.cos)
    self._compareBoth(
        y,
        np.vectorize(self._replace_domain_error_with_inf(math.lgamma)),
        tf.lgamma)
    self._compareBoth(x, np.vectorize(math.erf), tf.erf)
    self._compareBoth(x, np.vectorize(math.erfc), tf.erfc)
    self._compareBoth(x, np.arctan, tf.atan)
    self._compareBoth(k, np.arcsin, tf.asin)
    self._compareBoth(k, np.arccos, tf.acos)
    self._compareBoth(k, np.tan, tf.tan)

    self._compareBothSparse(x, np.abs, tf.abs)
    self._compareBothSparse(x, np.negative, tf.neg)
    self._compareBothSparse(x, np.square, tf.square)
    self._compareBothSparse(z, np.sqrt, tf.sqrt, tol=1e-3)
    self._compareBothSparse(x, np.tanh, tf.tanh)
    self._compareBothSparse(y, np.sign, tf.sign)
    self._compareBothSparse(x, np.vectorize(math.erf), tf.erf) 

def testHalfBasic(self):
    x = np.arange(-3, 3).reshape(1, 3, 2).astype(np.float16)
    y = (x + .5).astype(np.float16)    # no zero
    z = (x + 15.5).astype(np.float16)  # all positive
    self._compareBoth(x, np.abs, tf.abs)
    self._compareBoth(x, np.abs, _ABS)
    self._compareBoth(x, np.negative, tf.neg)
    self._compareBoth(x, np.negative, _NEG)
    self._compareBoth(y, self._inv, tf.inv)
    self._compareBoth(x, np.square, tf.square)
    self._compareBoth(z, np.sqrt, tf.sqrt)
    self._compareBoth(z, self._rsqrt, tf.rsqrt)
    self._compareBoth(x, np.exp, tf.exp)
    self._compareBoth(z, np.log, tf.log)
    self._compareBoth(z, np.log1p, tf.log1p)
    self._compareBoth(x, np.tanh, tf.tanh)
    self._compareBoth(x, self._sigmoid, tf.sigmoid)
    self._compareBoth(y, np.sign, tf.sign)
    self._compareBoth(x, np.sin, tf.sin)
    self._compareBoth(x, np.cos, tf.cos)
    self._compareBoth(
        y,
        np.vectorize(self._replace_domain_error_with_inf(math.lgamma)),
        tf.lgamma)
    self._compareBoth(x, np.vectorize(math.erf), tf.erf)
    self._compareBoth(x, np.vectorize(math.erfc), tf.erfc)

    self._compareBothSparse(x, np.abs, tf.abs)
    self._compareBothSparse(x, np.negative, tf.neg)
    self._compareBothSparse(x, np.square, tf.square)
    self._compareBothSparse(z, np.sqrt, tf.sqrt, tol=1e-3)
    self._compareBothSparse(x, np.tanh, tf.tanh)
    self._compareBothSparse(y, np.sign, tf.sign)
    self._compareBothSparse(x, np.vectorize(math.erf), tf.erf, tol=1e-3) 

def test_lgamma(self):
        tolerance = 14
        self.assertAlmostEqual(math.lgamma(0.5), 0.5 * math.log(math.pi), places=15)
        for i in range(1, 20):
            if i > 14:
                tolerance = 13
            self.assertAlmostEqual(math.log(math.factorial(i-1)), math.lgamma(i), places=tolerance)
        self.assertEqual(math.lgamma(float('inf')), float('inf'))
        self.assertEqual(math.lgamma(float('-inf')), float('inf'))
        self.assertTrue(math.isnan(math.lgamma(float('nan'))))
        for i in range(0, -1001, -1):
            self.assertRaises(ValueError, math.lgamma, i) 

def test_gamma_lgamma(self):
        tolerance = 13
        for x in itertools.count(0.001, 0.001):
            if x > 5.0:
                break
            self.assertAlmostEqual(math.lgamma(x), math.log(math.gamma(x)), places=tolerance)
            self.assertAlmostEqual(math.lgamma(x*x), math.log(math.gamma(x*x)), places=tolerance)
            self.assertAlmostEqual(math.lgamma(2.0**x), math.log(math.gamma(2.0**x)), places=tolerance)

            # Test negative values too, but not integers
            if x % 1.0 != 0.0:
                self.assertAlmostEqual(math.lgamma(-x), math.log(abs(math.gamma(-x))), places=tolerance)
                self.assertAlmostEqual(math.lgamma(-x*x), math.log(abs(math.gamma(-x*x))), places=tolerance)
                self.assertAlmostEqual(math.lgamma(-2.0**x), math.log(abs(math.gamma(-2.0**x))), places=tolerance) 

def _logfactorial(n):
  """Calculate natural logarithm of n!."""
  return math.lgamma(n + 1) 

def _gcf(a,x):
    "Continued fraction representation of Gamma. NumRec sect 6.1"
    ITMAX=100
    EPS=3.e-7
    FPMIN=1.e-30

    gln=lgamma(a)
    b=x+1.-a
    c=1./FPMIN
    d=1./b
    h=d
    for i in range(1,ITMAX+1):
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if abs(d) < FPMIN: d=FPMIN
        c=b+an/c
        if abs(c) < FPMIN: c=FPMIN
        d=1./d
        delt=d*c
        h=h*delt
        if abs(delt-1.) < EPS: break
    else:
        print('a too large, ITMAX too small in gcf')
    gammcf=np.exp(-x+a*np.log(x)-gln)*h
    return gammcf,gln 

def log_gampoiss(k,alpha,beta):
    import math
    k = int(k)
    return math.lgamma(k+alpha)+alpha*np.log(beta)-math.lgamma(alpha)-math.lgamma(k+1)-(alpha+k)*np.log(1+beta) 

def lgamma(x=('FloatPin', 0.0), Result=(REF, ('BoolPin', False))):
        '''Return the natural logarithm of the absolute value of the Gamma function at `x`.'''
        try:
            Result(True)
            return math.lgamma(x)
        except:
            Result(False)
            return -1 

def log_gamma(x):
    return math.lgamma(x) 

def logfactorial(n):

    return lgamma(n + 1) 

def multinomial_distance(p: np.ndarray, q: np.ndarray) -> float:
	N = p.shape[0]
	p_sum = p.sum()
	q_sum = q.sum()
	x = lgamma(N) + lgamma(p_sum + q_sum + N) - lgamma(p_sum + N) - lgamma(q_sum + N)
	for k in range(N):
		x += lgamma(p[k] + 1) + lgamma(q[k] + 1) - lgamma(1) - lgamma(p[k] + q[k] + 1)
	x = np.exp(x)
	return 1 - 1 / (1 + x)


# https://maciejkula.github.io/2015/02/22/incremental-construction-of-sparse-matrices/ 

def LogBetaBinomPMF(k, n, alpha, beta):
  lp = math.lgamma(n + 1)
  lp -= math.lgamma(k + 1) + math.lgamma(n - k + 1)
  lp += math.lgamma(k + alpha) + math.lgamma(n - k + beta)
  lp -= math.lgamma(n + alpha + beta)
  lp += math.lgamma(alpha + beta)
  lp -= math.lgamma(alpha) + math.lgamma(beta)
  return lp 

def testFloatBasic(self):
    x = np.arange(-3, 3).reshape(1, 3, 2).astype(np.float32)
    y = (x + .5).astype(np.float32)     # no zero
    z = (x + 15.5).astype(np.float32)   # all positive
    k = np.arange(-0.90, 0.90, 0.25).astype(np.float32) # between -1 and 1

    self._compareBoth(x, np.abs, tf.abs)
    self._compareBoth(x, np.abs, _ABS)
    self._compareBoth(x, np.negative, tf.neg)
    self._compareBoth(x, np.negative, _NEG)
    self._compareBoth(y, self._inv, tf.inv)
    self._compareBoth(x, np.square, tf.square)
    self._compareBoth(z, np.sqrt, tf.sqrt)
    self._compareBoth(z, self._rsqrt, tf.rsqrt)
    self._compareBoth(x, np.exp, tf.exp)
    self._compareBoth(z, np.log, tf.log)
    self._compareBoth(z, np.log1p, tf.log1p)
    self._compareBoth(x, np.tanh, tf.tanh)
    self._compareBoth(x, self._sigmoid, tf.sigmoid)
    self._compareBoth(y, np.sign, tf.sign)
    self._compareBoth(x, np.sin, tf.sin)
    self._compareBoth(x, np.cos, tf.cos)
    self._compareBoth(k, np.arcsin, tf.asin)
    self._compareBoth(k, np.arccos, tf.acos)
    self._compareBoth(x, np.arctan, tf.atan)
    self._compareBoth(x, np.tan, tf.tan)
    self._compareBoth(
        y,
        np.vectorize(self._replace_domain_error_with_inf(math.lgamma)),
        tf.lgamma)
    self._compareBoth(x, np.vectorize(math.erf), tf.erf)
    self._compareBoth(x, np.vectorize(math.erfc), tf.erfc)

    self._compareBothSparse(x, np.abs, tf.abs)
    self._compareBothSparse(x, np.negative, tf.neg)
    self._compareBothSparse(x, np.square, tf.square)
    self._compareBothSparse(z, np.sqrt, tf.sqrt, tol=1e-3)
    self._compareBothSparse(x, np.tanh, tf.tanh)
    self._compareBothSparse(y, np.sign, tf.sign)
    self._compareBothSparse(x, np.vectorize(math.erf), tf.erf) 

def test_pure_python_math_module(self):
        vals = [1, -.5, 1.5, 0, 0.0, -2, -2.2, .2]

        # not being tested: math.asinh, math.atanh, math.lgamma, math.erfc, math.acos
        def f():
            functions = [
                math.sqrt, math.cos, math.sin, math.tan, math.asin, math.atan,
                math.acosh, math.cosh, math.sinh, math.tanh, math.ceil,
                math.erf, math.exp, math.expm1, math.factorial, math.floor,
                math.log, math.log10, math.log1p
            ]
            tr = []
            for idx1 in range(len(vals)):
                v1 = vals[idx1]
                for funIdx in range(len(functions)):
                    function = functions[funIdx]
                    try:
                        tr = tr + [function(v1)]
                    except ValueError as ex:
                        pass

            return tr

        r1 = self.evaluateWithExecutor(f)
        r2 = f()
        self.assertGreater(len(r1), 100)
        self.assertTrue(numpy.allclose(r1, r2, 1e-6)) 

def __call__(self, val):
        if val <= -1 or val == 0:
            raise ValueError("math domain error")

        return __inline_fora(
            """fun(@unnamed_args:(val), *args) {
                   PyFloat(math.lgamma(val.@m))
                   }"""
            )(val) 

def test_lgamma(self):
        tolerance = 14
        self.assertAlmostEqual(math.lgamma(0.5), 0.5 * math.log(math.pi), places=15)
        for i in xrange(1, 20):
            if i > 14:
                tolerance = 13
            self.assertAlmostEqual(math.log(math.factorial(i-1)), math.lgamma(i), places=tolerance)
        self.assertEqual(math.lgamma(float('inf')), float('inf'))
        self.assertEqual(math.lgamma(float('-inf')), float('inf'))
        self.assertTrue(math.isnan(math.lgamma(float('nan'))))
        for i in xrange(0, -1001, -1):
            self.assertRaises(ValueError, math.lgamma, i) 

def lbeta(m, n):
    """Return log(Beta(m,n))"""
    return math.lgamma(m) + math.lgamma(n) - math.lgamma(m + n) 

def calc_log_Z(r, s, nu):
        return (
            ((nu + 1.) / 2.) * LOG2
            + .5 * LOGPI
            - .5 * log(r)
            - (nu/2.) * log(s)
            + lgamma(nu/2.)) 

def log_nCk(n, k):
    """log(choose(n,k)) with overflow protection."""
    if n == 0 or k == 0 or n == k:
        return 0
    return log(n) + lgamma(n) - log(k) - lgamma(k) - log(n-k) - lgamma(n-k)