Python math.lgamma方法代码示例

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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) 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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)) 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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) 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
def logp_crp(N, Nk, alpha):
    """Returns the log normalized P(N,K|alpha), where N is the number of
    customers and K is the number of tables.
    http://gershmanlab.webfactional.com/pubs/GershmanBlei12.pdf#page=4 (eq 8)
    """
    return len(Nk)*log(alpha) + np.sum(lgamma(c) for c in Nk) \
        + lgamma(alpha) - lgamma(N+alpha) 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
def logp_crp_unorm(N, K, alpha):
    """Returns the log unnormalized P(N,K|alpha), where N is the number of
    customers and K is the number of tables. Use for effeciency to avoid
    computing terms that are not a function of alpha.
    """
    return K*log(alpha) + lgamma(alpha) - lgamma(N+alpha) 

# 需要导入模块: import math [as 别名]
# 或者: from math import lgamma [as 别名]
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)