Python math 模块，lgamma() 实例源码

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)

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)

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)

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_gamma(x):
    return math.lgamma(x)

def log_betabin(k, n, alpha, beta):
    try:
        c = math.lgamma(alpha+beta) - math.lgamma(alpha) - math.lgamma(beta)
    except:
        print('alpha = {}, beta = {}'.format(alpha, beta))
    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(n-k+beta) - math.lgamma(n+alpha+beta) + c)

def log_betabin(k, n, alpha, beta):
    try:
        c = math.lgamma(alpha+beta) - math.lgamma(alpha) - math.lgamma(beta)
    except:
        print('alpha = {}, beta = {}'.format(alpha, beta))
    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(n-k+beta) - math.lgamma(n+alpha+beta) + c)

def kl_Beta(alpha, beta, alpha_0, beta_0):
    return tf.reduce_sum(math.lgamma(alpha_0) + math.lgamma(beta_0) - math.lgamma(alpha_0+beta_0)
                + tf.lgamma(alpha + beta) - tf.lgamma(alpha) - tf.lgamma(beta)
                + (alpha - alpha_0) * tf.digamma(alpha) + (beta - beta_0) * tf.digamma(beta)
                - (alpha + beta - alpha_0 - beta_0) * tf.digamma(alpha + beta)
        )

def nct(tval, delta, df):
  #This is the non-central t-distruction function
  #This is a direct translation from visual-basic to Python of the nct implementation by Geoff Cumming for ESCI

  #Noncentral t function, after Excel worksheet calculation by Geoff Robinson
  #The three parameters are:
  #   tval -- our t value
  #   delta -- noncentrality parameter
  #   df -- degrees of freedom

  #Dim dfmo As Integer     'for df-1
  #Dim tosdf As Double     'for tval over square root of df
  #Dim sep As Double       'for separation of points
  #Dim sepa As Double      'for temp calc of points
  #Dim consta As Double    'for constant
  #Dim ctr As Integer      'loop counter

  dfmo = df - 1.0
  tosdf = tval / math.sqrt(df)
  sep = (math.sqrt(df) + 7) / 100
  consta = math.exp( (2 - df) * 0.5 * math.log1p(1) - math.lgamma(df/2) ) * sep /3

  #now do first term in cross product summation, with df=0 special
  if dfmo > 0:
    nctvalue = stdnormdist(0 * tosdf - delta) * 0 ** dfmo * math.exp(-0.5* 0 * 0)
  else:  
    nctvalue = stdnormdist(0 * tosdf - delta) * math.exp(-0.5 * 0 * 0)

  #now add in second term, with multiplier 4
  nctvalue = nctvalue + 4 * stdnormdist(sep*tosdf - delta) * sep ** dfmo * math.exp(-0.5*sep*sep)    

  #now loop 49 times to add 98 terms
  for ctr in range(1,49):
      sepa = 2 * ctr * sep
      nctvalue = nctvalue + 2 * stdnormdist(sepa * tosdf - delta) * sepa ** dfmo * math.exp(-0.5 * sepa * sepa)
      sepa = sepa + sep
      nctvalue = nctvalue + 4 * stdnormdist(sepa * tosdf - delta) * sepa ** dfmo * math.exp(-0.5 * sepa * sepa)

  #add in last term
  sepa = sepa + sep
  nctvalue = nctvalue + stdnormdist(sepa * tosdf - delta) * sepa ** dfmo * math.exp(-0.5 * sepa * sepa)

  #multiply by the constant and we've finished
  nctvalue = nctvalue * consta

  return nctvalue