Python special.gammaln方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def test_value(self):
        with self.session(use_gpu=True):
            def _test_value(given, temperature, logits):
                given = np.array(given, np.float32)
                logits = np.array(logits, np.float32)
                n = logits.shape[-1]

                t = temperature
                target_log_p = gammaln(n) + (n - 1) * np.log(t) + \
                    (logits - t * given).sum(axis=-1) - \
                    n * np.log(np.exp(logits - t * given).sum(axis=-1))

                con = ExpConcrete(temperature, logits=logits)
                log_p = con.log_prob(given)
                self.assertAllClose(log_p.eval(), target_log_p)
                p = con.prob(given)
                self.assertAllClose(p.eval(), np.exp(target_log_p))

            _test_value([np.log(0.25), np.log(0.25), np.log(0.5)],
                        0.1,
                        [1., 1., 1.2])
            _test_value([[np.log(0.25), np.log(0.25), np.log(0.5)],
                        [np.log(0.1), np.log(0.5), np.log(0.4)]],
                        0.5,
                        [[1., 1., 1.], [.5, .5, .4]]) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def testGammaln(self):
        raw = np.random.rand(10, 8, 5)
        t = tensor(raw, chunk_size=3)

        r = gammaln(t)
        expect = scipy_gammaln(raw)

        self.assertEqual(r.shape, raw.shape)
        self.assertEqual(r.dtype, expect.dtype)

        r = r.tiles()
        t = get_tiled(t)

        self.assertEqual(r.nsplits, t.nsplits)
        for c in r.chunks:
            self.assertIsInstance(c.op, TensorGammaln)
            self.assertEqual(c.index, c.inputs[0].index)
            self.assertEqual(c.shape, c.inputs[0].shape) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def _pdf(self, x, df, nc):
        # nct.pdf(x, df, nc) =
        #                               df**(df/2) * gamma(df+1)
        #                ----------------------------------------------------
        #                2**df*exp(nc**2/2) * (df+x**2)**(df/2) * gamma(df/2)
        n = df*1.0
        nc = nc*1.0
        x2 = x*x
        ncx2 = nc*nc*x2
        fac1 = n + x2
        trm1 = n/2.*np.log(n) + sc.gammaln(n+1)
        trm1 -= n*np.log(2)+nc*nc/2.+(n/2.)*np.log(fac1)+sc.gammaln(n/2.)
        Px = np.exp(trm1)
        valF = ncx2 / (2*fac1)
        trm1 = np.sqrt(2)*nc*x*sc.hyp1f1(n/2+1, 1.5, valF)
        trm1 /= np.asarray(fac1*sc.gamma((n+1)/2))
        trm2 = sc.hyp1f1((n+1)/2, 0.5, valF)
        trm2 /= np.asarray(np.sqrt(fac1)*sc.gamma(n/2+1))
        Px *= trm1+trm2
        return Px 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def _lnB(alpha):
    r"""
    Internal helper function to compute the log of the useful quotient

    .. math::

        B(\alpha) = \frac{\prod_{i=1}{K}\Gamma(\alpha_i)}{\Gamma\left(\sum_{i=1}^{K}\alpha_i\right)}

    Parameters
    ----------
    %(_dirichlet_doc_default_callparams)s

    Returns
    -------
    B : scalar
        Helper quotient, internal use only

    """
    return np.sum(gammaln(alpha)) - gammaln(np.sum(alpha)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def lpol_fima(d, n=20):
    """MA representation of fractional integration

    .. math:: (1-L)^{-d} for |d|<0.5  or |d|<1 (?)

    Parameters
    ----------
    d : float
        fractional power
    n : int
        number of terms to calculate, including lag zero

    Returns
    -------
    ma : array
        coefficients of lag polynomial

    """
    # hide import inside function until we use this heavily
    from scipy.special import gammaln
    j = np.arange(n)
    return np.exp(gammaln(d + j) - gammaln(j + 1) - gammaln(d))


# moved from sandbox.tsa.try_fi 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def matrix_elem(n, r, m):
    """Matrix element corresponding to squeezed density matrix[n, m]"""
    eps = 1e-10

    if n % 2 != m % 2:
        return 0.0

    if r == 0.0:
        return np.complex(n == m)  # delta function

    k = np.arange(m % 2, min([m, n]) + 1, 2)
    res = np.sum(
        (-1) ** ((n - k) / 2)
        * np.exp(
            (lg(m + 1) + lg(n + 1)) / 2
            - lg(k + 1)
            - lg((m - k) / 2 + 1)
            - lg((n - k) / 2 + 1)
        )
        * (np.sinh(r) / 2 + eps) ** ((n + m - 2 * k) / 2)
        / (np.cosh(r) ** ((n + m + 1) / 2))
    )
    return res 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def lda_e_step(doc_word_ids, doc_word_counts, alpha, beta, max_iter=100):
    gamma = np.ones(len(alpha))
    expElogtheta = np.exp(dirichlet_expectation(gamma))
    betad = beta[:, doc_word_ids]
    phinorm = np.dot(expElogtheta, betad) + 1e-100
    counts = np.array(doc_word_counts)
    for _ in xrange(max_iter):
        lastgamma = gamma

        gamma = alpha + expElogtheta * np.dot(counts / phinorm, betad.T)
        Elogtheta = dirichlet_expectation(gamma)
        expElogtheta = np.exp(Elogtheta)
        phinorm = np.dot(expElogtheta, betad) + 1e-100
        meanchange = np.mean(abs(gamma - lastgamma))
        if (meanchange < meanchangethresh):
            break

    likelihood = np.sum(counts * np.log(phinorm))
    likelihood += np.sum((alpha - gamma) * Elogtheta)
    likelihood += np.sum(sp.gammaln(gamma) - sp.gammaln(alpha))
    likelihood += sp.gammaln(np.sum(alpha)) - sp.gammaln(np.sum(gamma))

    return (likelihood, gamma) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def _log_wishart_norm(degrees_of_freedom, log_det_precisions_chol, n_features):
    """Compute the log of the Wishart distribution normalization term.

    Parameters
    ----------
    degrees_of_freedom : array-like, shape (n_components,)
        The number of degrees of freedom on the covariance Wishart
        distributions.

    log_det_precision_chol : array-like, shape (n_components,)
         The determinant of the precision matrix for each component.

    n_features : int
        The number of features.

    Return
    ------
    log_wishart_norm : array-like, shape (n_components,)
        The log normalization of the Wishart distribution.
    """
    # To simplify the computation we have removed the np.log(np.pi) term
    return -(degrees_of_freedom * log_det_precisions_chol +
             degrees_of_freedom * n_features * .5 * math.log(2.) +
             np.sum(gammaln(.5 * (degrees_of_freedom -
                                  np.arange(n_features)[:, np.newaxis])), 0)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def test_log_wishart_norm():
    rng = np.random.RandomState(0)

    n_components, n_features = 5, 2
    degrees_of_freedom = np.abs(rng.rand(n_components)) + 1.
    log_det_precisions_chol = n_features * np.log(range(2, 2 + n_components))

    expected_norm = np.empty(5)
    for k, (degrees_of_freedom_k, log_det_k) in enumerate(
            zip(degrees_of_freedom, log_det_precisions_chol)):
        expected_norm[k] = -(
            degrees_of_freedom_k * (log_det_k + .5 * n_features * np.log(2.)) +
            np.sum(gammaln(.5 * (degrees_of_freedom_k -
                                 np.arange(0, n_features)[:, np.newaxis])), 0))
    predected_norm = _log_wishart_norm(degrees_of_freedom,
                                       log_det_precisions_chol, n_features)

    assert_almost_equal(expected_norm, predected_norm) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def meanfieldupdate_p(self, stepsize=1.0):
        """
        Update p given the network parameters and the current variational
        parameters of the weight distributions.
        :return:
        """
        logit_p = self.network.expected_log_p() - self.network.expected_log_notp()
        logit_p += self.network.kappa * self.network.expected_log_v() - gammaln(self.network.kappa)
        logit_p += gammaln(self.mf_kappa_1) - self.mf_kappa_1 * np.log(self.mf_v_1)
        logit_p += gammaln(self.kappa_0) - self.kappa_0 * np.log(self.nu_0)
        logit_p += self.mf_kappa_0 * np.log(self.mf_v_0) - gammaln(self.mf_kappa_0)

        # p_hat = logistic(logit_p)
        # self.mf_p = (1.0 - stepsize) * self.mf_p + stepsize * p_hat

        logit_p_hat = (1-stepsize) * logit(self.mf_p) + \
                       stepsize * logit_p
        # self.mf_p = logistic(logit_p_hat)
        self.mf_p = np.clip(logistic(logit_p_hat), 1e-8, 1-1e-8) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def _gamma(N):
    # Note: this is closely approximated by (1 - 0.75 / N) for large N
    return np.sqrt(2 / N) * np.exp(gammaln(N / 2) - gammaln((N - 1) / 2)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def _log_gamma(N):
    return 0.5 * np.log(2 / N) + gammaln(N / 2) - gammaln((N - 1) / 2) 

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

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def incorporate(self, rowid, observation, inputs=None):
        DistributionGpm.incorporate(self, rowid, observation, inputs)
        x = observation[self.outputs[0]]
        if not (x % 1 == 0 and x >= 0):
            raise ValueError('Invalid Poisson: %s' % str(x))
        self.N += 1
        self.sum_x += x
        self.sum_log_fact_x += gammaln(x+1)
        self.data[rowid] = x 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import gammaln [as 别名]
def unincorporate(self, rowid):
        x = self.data.pop(rowid)
        self.N -= 1
        self.sum_x -= x
        self.sum_log_fact_x -= gammaln(x+1)