Python special.betaln方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def AS_betabinom_loglike(logps, sigma, AS1, AS2, hetp, error):
    if sigma >= 1.0 or sigma <= 0.0:
        return -99999999999.0

    a = math.exp(logps[0] + math.log(1/sigma**2 - 1))
    b = math.exp(logps[1] + math.log(1/sigma**2 - 1))

    part1 = 0
    part1 += betaln(AS1 + a, AS2 + b)
    part1 -= betaln(a, b)

    if hetp == 1:
        return part1

    e1 = math.log(error) * AS1 + math.log(1 - error) * AS2
    e2 = math.log(error) * AS2 + math.log(1 - error) * AS1
    if hetp == 0:
        return addlogs(e1, e2)

    return addlogs(math.log(hetp)+part1, math.log(1-hetp) + addlogs(e1, e2)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def AS_betabinom_loglike(logps, sigma, AS1, AS2, hetp, error):
    """Given parameter, returns log likelihood of allele-specific
    part of test. Note that some parts of equation have been
    canceled out"""
    a = math.exp(logps[0] + math.log(1/sigma**2 - 1))
    b = math.exp(logps[1] + math.log(1/sigma**2 - 1))

    part1 = 0
    part1 += betaln(AS1 + a, AS2 + b)
    part1 -= betaln(a, b)

    if hetp==1:
        return part1

    e1 = math.log(error) * AS1 + math.log(1 - error) * AS2
    e2 = math.log(error) * AS2 + math.log(1 - error) * AS1
    if hetp == 0:
        return addlogs(e1, e2)

    return addlogs(math.log(hetp)+part1, math.log(1-hetp) + addlogs(e1, e2)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def calc_log_Z(a, b):
        return betaln(a, b) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def calc_logpdf_marginal(N, x_sum, alpha, beta):
        return betaln(x_sum + alpha, N - x_sum + beta) - betaln(alpha, beta) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def _logpdf(self, x, a, b):
        lPx = sc.xlog1py(b - 1.0, -x) + sc.xlogy(a - 1.0, x)
        lPx -= sc.betaln(a, b)
        return lPx 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def _logpmf(self, k, M, n, N):
        tot, good = M, n
        bad = tot - good
        return betaln(good+1, 1) + betaln(bad+1,1) + betaln(tot-N+1, N+1)\
            - betaln(k+1, good-k+1) - betaln(N-k+1,bad-N+k+1)\
            - betaln(tot+1, 1) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def test_betaln(self):
        assert_equal(cephes.betaln(1,1),0.0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def test_beta():
    np.random.seed(1234)

    b = np.r_[np.logspace(-200, 200, 4),
              np.logspace(-10, 10, 4),
              np.logspace(-1, 1, 4),
              np.arange(-10, 11, 1),
              np.arange(-10, 11, 1) + 0.5,
              -1, -2.3, -3, -100.3, -10003.4]
    a = b

    ab = np.array(np.broadcast_arrays(a[:,None], b[None,:])).reshape(2, -1).T

    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 400

        assert_func_equal(sc.beta,
                          lambda a, b: float(mpmath.beta(a, b)),
                          ab,
                          vectorized=False,
                          rtol=1e-10,
                          ignore_inf_sign=True)

        assert_func_equal(
            sc.betaln,
            lambda a, b: float(mpmath.log(abs(mpmath.beta(a, b)))),
            ab,
            vectorized=False,
            rtol=1e-10)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


#------------------------------------------------------------------------------
# Machinery for systematic tests
#------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def test_betaln(self):
        assert_equal(cephes.betaln(1,1),0.0)
        assert_allclose(cephes.betaln(-100.3, 1e-200), cephes.gammaln(1e-200))
        assert_allclose(cephes.betaln(0.0342, 170), 3.1811881124242447,
                        rtol=1e-14, atol=0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def test_beta():
    np.random.seed(1234)

    b = np.r_[np.logspace(-200, 200, 4),
              np.logspace(-10, 10, 4),
              np.logspace(-1, 1, 4),
              np.arange(-10, 11, 1),
              np.arange(-10, 11, 1) + 0.5,
              -1, -2.3, -3, -100.3, -10003.4]
    a = b

    ab = np.array(np.broadcast_arrays(a[:,None], b[None,:])).reshape(2, -1).T

    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 400

        assert_func_equal(sc.beta,
                          lambda a, b: float(mpmath.beta(a, b)),
                          ab,
                          vectorized=False,
                          rtol=1e-10,
                          ignore_inf_sign=True)

        assert_func_equal(
            sc.betaln,
            lambda a, b: float(mpmath.log(abs(mpmath.beta(a, b)))),
            ab,
            vectorized=False,
            rtol=1e-10)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


# ------------------------------------------------------------------------------
# loggamma
# ------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def testBetaBetaKL(self):
    with self.test_session() as sess:
      for shape in [(10,), (4,5)]:
        a1 = 6.0*np.random.random(size=shape) + 1e-4
        b1 = 6.0*np.random.random(size=shape) + 1e-4 
        a2 = 6.0*np.random.random(size=shape) + 1e-4
        b2 = 6.0*np.random.random(size=shape) + 1e-4 
        # Take inverse softplus of values to test BetaWithSoftplusAB
        a1_sp = np.log(np.exp(a1) - 1.0)
        b1_sp = np.log(np.exp(b1) - 1.0)
        a2_sp = np.log(np.exp(a2) - 1.0)
        b2_sp = np.log(np.exp(b2) - 1.0)

        d1 = tf.contrib.distributions.Beta(a=a1, b=b1)
        d2 = tf.contrib.distributions.Beta(a=a2, b=b2)
        d1_sp = tf.contrib.distributions.BetaWithSoftplusAB(a=a1_sp, b=b1_sp)
        d2_sp = tf.contrib.distributions.BetaWithSoftplusAB(a=a2_sp, b=b2_sp)

        kl_expected = (special.betaln(a2, b2) - special.betaln(a1, b1)
                     + (a1 - a2)*special.digamma(a1)
                     + (b1 - b2)*special.digamma(b1)
                     + (a2 - a1 + b2 - b1)*special.digamma(a1 + b1))

        for dist1 in [d1, d1_sp]:
          for dist2 in [d2, d2_sp]:
            kl = tf.contrib.distributions.kl(dist1, dist2)
            kl_val = sess.run(kl)
            self.assertEqual(kl.get_shape(), shape)
            self.assertAllClose(kl_val, kl_expected)
        
        # Make sure KL(d1||d1) is 0
        kl_same = sess.run(tf.contrib.distributions.kl(d1, d1))
        self.assertAllClose(kl_same, np.zeros_like(kl_expected)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def lbeta_asymp(a, b):
    if b > a:
        a, b = b, a

    if a<1e6:
        return betaln(a, b)

    l = gammaln(b)

    l -= b*math.log(a)
    l += b*(1-b)/(2*a)
    l += b*(1-b)*(1-2*b)/(12*a*a)
    l += -((b*(1-b))**2)/(12*a**3)

    return l 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def BNB_loglike(k, mean, n, sigma):
    #n=min(n,10000)
    #Put variables in beta-NB form (n,a,b)
    mean = max(mean, 0.00001)
    p = np.float64(n)/(n+mean)
    logps = [math.log(n) - math.log(n + mean),
             math.log(mean) - math.log(n + mean)]

    if sigma < 0.00001:
        loglike = -betaln(n, k+1)-math.log(n+k)+n*logps[0]+k*logps[1]
        return loglike

    sigma = (1/sigma)**2

    a = p * sigma + 1
    b = (1-p) * sigma

    loglike = 0

    #Rising Pochhammer = gamma(k+n)/gamma(n)
    #for j in range(k):
    #    loglike += math.log(j+n)
    if k>0:
        loglike = -lbeta_asymp(n, k) - math.log(k)
        #loglike=scipy.special.gammaln(k+n)-scipy.special.gammaln(n)
    else:
        loglike=0

    #Add log(beta(a+n,b+k))
    loglike += lbeta_asymp(a+n, b+k)

    #Subtract log(beta(a,b))
    loglike -= lbeta_asymp(a, b)

    return loglike 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def betaln_asym(a, b):
    if b > a:
        a, b = b, a

    if a < 1e6:
        return betaln(a, b)

    l=gammaln(b)

    l -= b*math.log(a)
    l += b*(1-b)/(2*a)
    l += b*(1-b)*(1-2*b)/(12*a*a)
    l += -((b*(1-b))**2)/(12*a**3)

    return l 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import betaln [as 别名]
def dominance(variations, control=None, sample_size=SAMPLE_SIZE):
    """Calculate P(A > B).

    Uses a modified Evan Miller closed formula if prior parameters are integers
    http://www.evanmiller.org/bayesian-ab-testing.html

    Uses scipy's MCMC otherwise

    TODO: The modified formula for informative prior has to be proved correct
    """
    values = OrderedDict()
    a, others = _split(variations, control)

    def is_integer(x):
        try:
            return x.is_integer()
        except:
            return int(x) == x

    for label, b in others.items():

        # If prior parameters are integers, use modified Evan Miller formula:
        if is_integer(a.prior_alpha) and is_integer(a.prior_beta) \
           and is_integer(b.prior_alpha) and is_integer(b.prior_beta):

            total = 0
            for i in range(b.alpha-b.prior_alpha):
                total += np.exp(spc.betaln(a.alpha+i, b.beta + a.beta) -
                                np.log(b.beta+i) - spc.betaln(b.prior_alpha+i,
                                b.beta) - spc.betaln(a.alpha, a.beta))
            values[label] = total

        # Use MCMC otherwise
        else:
            a_samples = a.posterior.rvs(sample_size)
            b_samples = b.posterior.rvs(sample_size)
            values[label] = np.mean(b_samples > a_samples)

    return values