Python math.lgamma() Examples

The following are code examples for showing how to use math.lgamma(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: pydistinct   Author: chanedwin   File: utils.py    MIT License 6 votes vote down vote up
def memoized_gamma(x, memo_dict):
    """
    standard implementation of the gamma function

    :param x:value to evaluate gamma function at
    :type x: float
    :param memo_dict: memoized dictionary for precomputed gamma values
    :type memo_dict: dict
    :return: value of gamma function evaluated at x, and memoized results
    :rtype: int, dict
    """
    x = int(x)
    if x in memo_dict:
        return memo_dict[x], memo_dict
    else:
        result = math.lgamma(x)
        memo_dict[x] = result
        return result, memo_dict 
Example 2
Project: ad_examples   Author: shubhomoydas   File: bayesian_ruleset.py    MIT License 6 votes vote down vote up
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 
Example 3
Project: quadpy   Author: nschloe   File: _helpers.py    MIT License 6 votes vote down vote up
def integrate_monomial_over_unit_nsphere(alpha, symbolic=False):
    """
    Gerald B. Folland,
    How to Integrate a Polynomial over a Sphere,
    The American Mathematical Monthly,
    Vol. 108, No. 5 (May, 2001), pp. 446-448,
    <https://doi.org/10.2307/2695802>.
    """
    if any(a % 2 == 1 for a in alpha):
        return 0

    if symbolic:
        return 2 * (
            prod([gamma(Rational(a + 1, 2)) for a in alpha])
            / gamma(sum([Rational(a + 1, 2) for a in alpha]))
        )

    # Use lgamma since other with ordinary gamma, numerator and denominator
    # might overflow.
    return 2 * math.exp(
        math.fsum([math.lgamma(0.5 * (a + 1)) for a in alpha])
        - math.lgamma(math.fsum([0.5 * (a + 1) for a in alpha]))
    ) 
Example 4
Project: molecular-design-toolkit   Author: Autodesk   File: utils.py    Apache License 2.0 6 votes vote down vote up
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 
Example 5
Project: titus2   Author: animator   File: dist.py    Apache License 2.0 6 votes vote down vote up
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)) 
Example 6
Project: titus2   Author: animator   File: dist.py    Apache License 2.0 6 votes vote down vote up
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 
Example 7
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 6 votes vote down vote up
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 
Example 8
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 6 votes vote down vote up
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 
Example 9
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 6 votes vote down vote up
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 
Example 10
Project: BOA   Author: wangtongada   File: BOAmodel.py    MIT License 6 votes vote down vote up
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 
Example 11
Project: ironpython2   Author: IronLanguages   File: test_math.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 12
Project: ironpython2   Author: IronLanguages   File: test_math.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 13
Project: nucleus   Author: google   File: genomics_math.py    Apache License 2.0 5 votes vote down vote up
def log10_binomial(k, n, p):
  """Calculates numerically-stable value of log10(binomial(k, n, p)).

  Returns the log10 of the binomial density for k successes in n trials where
  each success has a probability of occurring of p.

  In real-space, we would calculate:

     result = (n choose k) * (1-p)^(n-k) * p^k

  This function computes the log10 of result, which is:

     log10(result) = log10(n choose k) + (n-k) * log10(1-p) + k * log10(p)

  This is equivalent to invoking the R function:
    dbinom(x=k, size=n, prob=p, log=TRUE)

  See https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Binomial.html
  for more details on the binomial.

  Args:
    k: int >= 0. Number of successes.
    n: int >= k. Number of trials.
    p: 0.0 <= float <= 1.0. Probability of success.

  Returns:
    log10 probability of seeing k successes in n trials with p.
  """
  r = math.lgamma(n + 1) - (math.lgamma(k + 1) + math.lgamma(n - k + 1))
  if k > 0:
    r += k * math.log(p)
  if n > k:
    r += (n-k) * math.log1p(-p)
  return r / LOG_E_OF_10 
Example 14
Project: gpo-smc-abc   Author: compops   File: models_dists.py    GNU General Public License v3.0 5 votes vote down vote up
def multiTLogPDF(x, mu, Sigma, nu, p):

    part1 = math.lgamma(0.5 * (p + nu))
    part2 = - math.lgamma(0.5 * nu) - 0.5 * p * \
        np.log(nu) - 0.5 * p * np.log(np.pi)
    part3 = - 0.5 * np.log(np.linalg.det(Sigma))
    part4 = - 0.5 * (nu + p) * np.log(1.0 + nu**(-1) *
                                      np.dot(np.dot((x - mu), np.linalg.inv(Sigma)), (x - mu)))

    return part1 + part2 + part3 + part4 
Example 15
Project: pynam   Author: hbp-unibi   File: entropy.py    GNU General Public License v3.0 5 votes vote down vote up
def lnncrr(x, y):
    """
    Returns the natural logarithm of the binomial coefficient for real
    arguments.
    """
    return (math.lgamma(x+1.0) - math.lgamma(y+1.0) - math.lgamma(x-y+1.0)) 
Example 16
Project: ufora   Author: ufora   File: pure_math.py    Apache License 2.0 5 votes vote down vote up
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([email protected]))
                   }"""
            )(val) 
Example 17
Project: ufora   Author: ufora   File: MathTestCases.py    Apache License 2.0 5 votes vote down vote up
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)) 
Example 18
Project: deep_image_model   Author: tobegit3hub   File: cwise_ops_test.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 19
Project: deep_image_model   Author: tobegit3hub   File: cwise_ops_test.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 20
Project: deep_image_model   Author: tobegit3hub   File: cwise_ops_test.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 21
Project: deep_image_model   Author: tobegit3hub   File: cwise_ops_test.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 22
Project: transform   Author: tensorflow   File: info_theory.py    Apache License 2.0 5 votes vote down vote up
def _logfactorial(n):
  """Calculate natural logarithm of n!."""
  return math.lgamma(n + 1) 
Example 23
Project: PyOpenDial   Author: KAIST-AILab   File: math_utils.py    MIT License 5 votes vote down vote up
def log_gamma(x):
        """
        Returns the log-gamma value using Lanczos approximation formula

        :param x: the point
        :return: the log-gamma value for the point
        """
        return math.lgamma(x) 
Example 24
Project: molecular-design-toolkit   Author: Autodesk   File: utils.py    Apache License 2.0 5 votes vote down vote up
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 
Example 25
Project: titus2   Author: animator   File: dist.py    Apache License 2.0 5 votes vote down vote up
def PDF(self, x):
        if (self.lamda == 0.0):
            if (x != 0.0):
                return 0.0
            else:
                return 1.0
        elif (x < 0.0):
            return 0.0
        else:
            return math.exp(x * math.log(self.lamda) - self.lamda - math.lgamma(x + 1.0)) 
Example 26
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 5 votes vote down vote up
def __call__(self, state, scope, pos, paramTypes, a):
        if math.isinf(a) and a > 0:
            return float("nan")
        try:
            if a <= 0: raise Exception
            else:
                return math.lgamma(a)
        except:
            raise PFARuntimeException("domain error", self.errcodeBase + 0, self.name, pos) 
Example 27
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 5 votes vote down vote up
def logBetaFunction(a, b):
    if (a <= 0.0) or (b <= 0.0):
        raise Exception
    else:
        return math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b)

# numerical recipes in C 
Example 28
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 5 votes vote down vote up
def gcf(a, x, itmax=200, eps=3.0e-7):
    """Continued fraction approximation of the incomplete gamma function."""
    gln = math.lgamma(a)
    gold = 0.0
    a0 = 1.0
    a1 = x
    b0 = 0.0
    b1 = 1.0
    fac = 1.0
    n = 1
    while n <= itmax:
        an = n
        ana = an - a
        a0 = (a1 + a0*ana)*fac
        b0 = (b1 + b0*ana)*fac
        anf = an*fac
        a1 = x*a0 + anf*a1
        b1 = x*b0 + anf*b1
        if (a1 != 0.0):
            fac = 1.0 / a1
            g = b1*fac
            if (abs((g-gold)/g) < eps):
                return (g*math.exp(-x+a*math.log(x)-gln), gln)
            gold = g
        n = n + 1
    raise Exception 
Example 29
Project: BOA   Author: wangtongada   File: BOAmodel.py    MIT License 5 votes vote down vote up
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) 
Example 30
Project: PyFlow   Author: wonderworks-software   File: MathLib.py    Apache License 2.0 5 votes vote down vote up
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 
Example 31
Project: m-stp   Author: MukeunKim   File: core.py    MIT License 5 votes vote down vote up
def lgamma(x):
    if isinstance(x, LineValue): lx = x.get_value()
    else: lx = x
    if lx == NAN: return LineValue(NAN)
    return LineValue(math.lgamma(lx))
# endregion 
Example 32
Project: cms   Author: broadinstitute   File: stats.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def log_choose(n, k) :
    # Return log(n choose k). Compute using lgamma(x + 1) = log(x!)
    if not (0 <= k <=n) :
        raise ValueError('%d is negative or more than %d' % (k, n))
    return lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1) 
Example 33
Project: Turing   Author: TuringApp   File: stats.py    MIT License 5 votes vote down vote up
def log_gamma(x):
    return math.lgamma(x) 
Example 34
Project: rtchange   Author: hyusuk   File: coding.py    MIT License 5 votes vote down vote up
def length(self, x):
        """Calculate SDNML code-length.

        Parameters
        ----------
        x : float
            Sample.

        Returns
        -------
        code_length : float
            Code length of the sample.
        """
        self._xs.pop(0)
        self._xs.append(x)
        xs = np.atleast_2d(self._xs).T
        self._update_stats(xs)
        code_length = 0
        if self._time > 1:
            code_length = (math.log(math.pi)/2) \
                          - math.log(zero2small(1-self._stats['d']))
            code_length -= math.lgamma((self._time-1)/2)
            code_length += math.lgamma((self._time)/2)
            code_length += (math.log(self._time-1) +
                            math.log(zero2small(self._stats['prev_tau'])))/2
            code_length += (self._time)/2 * (
                math.log(zero2small((self._time)*self._stats['tau'])) -
                math.log(zero2small((self._time-1)*self._stats['prev_tau']))
            )
        self._time += 1
        return code_length 
Example 35
Project: data_algebra   Author: WinVector   File: SQLite.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def prepare_connection(self, conn):
        # https://docs.python.org/3/library/sqlite3.html#sqlite3.Connection.create_function
        conn.create_function("is_bad", 1, _check_scalar_bad)
        # math fns
        conn.create_function("acos", 1, math.acos)
        conn.create_function("acosh", 1, math.acosh)
        conn.create_function("asin", 1, math.asin)
        conn.create_function("asinh", 1, math.asinh)
        conn.create_function("atan", 1, math.atan)
        conn.create_function("atanh", 1, math.atanh)
        conn.create_function("ceil", 1, math.ceil)
        conn.create_function("cos", 1, math.cos)
        conn.create_function("cosh", 1, math.cosh)
        conn.create_function("degrees", 1, math.degrees)
        conn.create_function("erf", 1, math.erf)
        conn.create_function("erfc", 1, math.erfc)
        conn.create_function("exp", 1, math.exp)
        conn.create_function("expm1", 1, math.expm1)
        conn.create_function("fabs", 1, math.fabs)
        conn.create_function("factorial", 1, math.factorial)
        conn.create_function("floor", 1, math.floor)
        conn.create_function("frexp", 1, math.frexp)
        conn.create_function("gamma", 1, math.gamma)
        conn.create_function("isfinite", 1, math.isfinite)
        conn.create_function("isinf", 1, math.isinf)
        conn.create_function("isnan", 1, math.isnan)
        conn.create_function("lgamma", 1, math.lgamma)
        conn.create_function("log", 1, math.log)
        conn.create_function("log10", 1, math.log10)
        conn.create_function("log1p", 1, math.log1p)
        conn.create_function("log2", 1, math.log2)
        conn.create_function("modf", 1, math.modf)
        conn.create_function("radians", 1, math.radians)
        conn.create_function("sin", 1, math.sin)
        conn.create_function("sinh", 1, math.sinh)
        conn.create_function("sqrt", 1, math.sqrt)
        conn.create_function("tan", 1, math.tan)
        conn.create_function("tanh", 1, math.tanh)
        conn.create_function("trunc", 1, math.trunc)
        conn.create_function("atan2", 2, math.atan2)
        conn.create_function("copysign", 2, math.copysign)
        conn.create_function("fmod", 2, math.fmod)
        conn.create_function("gcd", 2, math.gcd)
        conn.create_function("hypot", 2, math.hypot)
        conn.create_function("isclose", 2, math.isclose)
        conn.create_function("ldexp", 2, math.ldexp)
        conn.create_function("pow", 2, math.pow) 
Example 36
Project: data_algebra   Author: WinVector   File: SQLite.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def prepare_connection(self, conn):
        # https://docs.python.org/3/library/sqlite3.html#sqlite3.Connection.create_function
        conn.create_function("is_bad", 1, _check_scalar_bad)
        # math fns
        conn.create_function("acos", 1, math.acos)
        conn.create_function("acosh", 1, math.acosh)
        conn.create_function("asin", 1, math.asin)
        conn.create_function("asinh", 1, math.asinh)
        conn.create_function("atan", 1, math.atan)
        conn.create_function("atanh", 1, math.atanh)
        conn.create_function("ceil", 1, math.ceil)
        conn.create_function("cos", 1, math.cos)
        conn.create_function("cosh", 1, math.cosh)
        conn.create_function("degrees", 1, math.degrees)
        conn.create_function("erf", 1, math.erf)
        conn.create_function("erfc", 1, math.erfc)
        conn.create_function("exp", 1, math.exp)
        conn.create_function("expm1", 1, math.expm1)
        conn.create_function("fabs", 1, math.fabs)
        conn.create_function("factorial", 1, math.factorial)
        conn.create_function("floor", 1, math.floor)
        conn.create_function("frexp", 1, math.frexp)
        conn.create_function("gamma", 1, math.gamma)
        conn.create_function("isfinite", 1, math.isfinite)
        conn.create_function("isinf", 1, math.isinf)
        conn.create_function("isnan", 1, math.isnan)
        conn.create_function("lgamma", 1, math.lgamma)
        conn.create_function("log", 1, math.log)
        conn.create_function("log10", 1, math.log10)
        conn.create_function("log1p", 1, math.log1p)
        conn.create_function("log2", 1, math.log2)
        conn.create_function("modf", 1, math.modf)
        conn.create_function("radians", 1, math.radians)
        conn.create_function("sin", 1, math.sin)
        conn.create_function("sinh", 1, math.sinh)
        conn.create_function("sqrt", 1, math.sqrt)
        conn.create_function("tan", 1, math.tan)
        conn.create_function("tanh", 1, math.tanh)
        conn.create_function("trunc", 1, math.trunc)
        conn.create_function("atan2", 2, math.atan2)
        conn.create_function("copysign", 2, math.copysign)
        conn.create_function("fmod", 2, math.fmod)
        conn.create_function("gcd", 2, math.gcd)
        conn.create_function("hypot", 2, math.hypot)
        conn.create_function("isclose", 2, math.isclose)
        conn.create_function("ldexp", 2, math.ldexp)
        conn.create_function("pow", 2, math.pow) 
Example 37
Project: titus2   Author: animator   File: spec.py    Apache License 2.0 4 votes vote down vote up
def inverseIncompleteBetaFunction(p, a, b):
    a1 = a - 1.0
    b1 = b - 1.0
    ERROR = 1.0e-10

    if (p <= 0.0):
        return 0.0
    elif (p >= 1.0):
        return 1.0;
    elif (a >= 1.0) and (b >= 1.0):
        if (p < 0.5):
            pp = p
        else:
            pp = 1.0 - p;
        t = math.sqrt(-2.0*math.log(pp))
        x = (2.30753+t*0.27061)/(1.0+t*(0.99229+t*0.04481)) - t
        if (p < 0.5):
            x = -x
        al = ((x*x)-3.0)/6.0;
        h = 2.0/(1.0/(2.0*a-1.0)+1.0/(2.0*b-1.0));
        w = (x*math.sqrt(al+h)/h)-(1.0/(2.0*b-1)-1.0/(2.0*a-1.0))*(al+5.0/6.0-2.0/(3.0*h))
        x = a/(a+b*math.exp(2.0*w))
    else:
        lna = math.log(a/(a+b))
        lnb = math.log(b/(a+b))
        t = math.exp(a*lna)/a
        u = math.exp(b*lnb)/b
        w = t + u
        if (p < t/w):
            x = math.pow(a*w*p,1.0/a)
        else:
            x = 1.0 - math.pow(b*w*(1.0-p), 1.0/b)
    afac = -math.lgamma(a)-math.lgamma(b)+math.lgamma(a+b)
    j = 0
    for i in range(10):
        if (x == 0.0) or (x == 1.0):
            return x;
        err = incompleteBetaFunction(x,a,b) - p

        t = math.exp(a1*math.log(x)+b1*math.log(1.0-x) + afac)
        u = err/t
        t = u/(1.0-0.5*min(1.0,u*(a1/x - b1/(1.0-x))))
        x -= t
        if (x <= 0.0):
            x = 0.5*(x + t)
        if (x >= 1.0):
            x = 0.5*(x + t + 1.0)
        if (math.fabs(t) < ERROR*x) and (j > 0):
            break
    return x

# numerical recipes in C
# used by incompleteBeta and incompleteBetaComp