Python special.loggamma方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def KL_davidson(k, d):
    vmf_entropy = k * ive(d / 2, k) / ive((d / 2) - 1, k) + \
                  (d / 2 - 1) * np.log(k) \
                  - (d / 2) * np.log(2 * np.pi) - np.log(iv(d / 2 - 1, k))

    hyu_ent = np.log(2) + (d / 2) * np.log(np.pi) - sp.loggamma(
        d / 2)

    kl = vmf_entropy + hyu_ent
    return kl
#
# first = k * bessel(d / 2, k) / bessel(d / 2 - 1, k)
# second = (d / 2 - 1) * torch.log(k) - torch.log(bessel(d / 2 - 1, k))
# const = torch.tensor(
#            np.log(3.1415926) * d / 2 + np.log(2) - sp.loggamma(d / 2).real - (d / 2) * np.log(2 * 3.1415926)).to(
#            devic

# for kappa in range(10, 150, 20):
#     for d in range(50, 150, 50):
#         print("Davidson:{}\t\tGuu:{}".format(KL_davidson(kappa, d), KL_guu(kappa, d))) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def test_loggamma_taylor_transition():
    # Make sure there isn't a big jump in accuracy when we move from
    # using the Taylor series to using the recurrence relation.

    r = LOGGAMMA_TAYLOR_RADIUS + np.array([-0.1, -0.01, 0, 0.01, 0.1])
    theta = np.linspace(0, 2*np.pi, 20)
    r, theta = np.meshgrid(r, theta)
    dz = r*np.exp(1j*theta)
    z = np.r_[1 + dz, 2 + dz].flatten()

    dataset = []
    for z0 in z:
        dataset.append((z0, complex(mpmath.loggamma(z0))))
    dataset = np.array(dataset)

    FuncData(sc.loggamma, dataset, 0, 1, rtol=5e-14).check() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def test_loggamma_taylor():
    # Test around the zeros at z = 1, 2.

    r = np.logspace(-16, np.log10(LOGGAMMA_TAYLOR_RADIUS), 10)
    theta = np.linspace(0, 2*np.pi, 20)
    r, theta = np.meshgrid(r, theta)
    dz = r*np.exp(1j*theta)
    z = np.r_[1 + dz, 2 + dz].flatten()

    dataset = []
    for z0 in z:
        dataset.append((z0, complex(mpmath.loggamma(z0))))
    dataset = np.array(dataset)

    FuncData(sc.loggamma, dataset, 0, 1, rtol=5e-14).check()


# ------------------------------------------------------------------------------
# rgamma
# ------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def spectral_density(self, k):  # noqa: D102
        k = np.array(k, dtype=np.double)
        # for nu > 20 we just use an approximation of the gaussian model
        if self.nu > 20.0:
            return (
                (self.len_scale / np.sqrt(np.pi)) ** self.dim
                * np.exp(-((k * self.len_scale) ** 2))
                * (
                    1
                    + (
                        ((k * self.len_scale) ** 2 - self.dim / 2.0) ** 2
                        - self.dim / 2.0
                    )
                    / self.nu
                )
            )
        return (self.len_scale / np.sqrt(np.pi)) ** self.dim * np.exp(
            -(self.nu + self.dim / 2.0)
            * np.log(1.0 + (k * self.len_scale) ** 2 / self.nu)
            + sps.loggamma(self.nu + self.dim / 2.0)
            - sps.loggamma(self.nu)
            - self.dim * np.log(np.sqrt(self.nu))
        ) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def _lazywhere(cond, arrays, f, fillvalue=None, f2=None):
    """
    np.where(cond, x, fillvalue) always evaluates x even where cond is False.
    This one only evaluates f(arr1[cond], arr2[cond], ...).
    For example,
    >>> a, b = np.array([1, 2, 3, 4]), np.array([5, 6, 7, 8])
    >>> def f(a, b):
        return a*b
    >>> _lazywhere(a > 2, (a, b), f, np.nan)
    array([ nan,  nan,  21.,  32.])
    Notice it assumes that all `arrays` are of the same shape, or can be
    broadcasted together.
    """
    if fillvalue is None:
        if f2 is None:
            raise ValueError("One of (fillvalue, f2) must be given.")
        else:
            fillvalue = np.nan
    else:
        if f2 is not None:
            raise ValueError("Only one of (fillvalue, f2) can be given.")

    arrays = np.broadcast_arrays(*arrays)
    temp = tuple(np.extract(cond, arr) for arr in arrays)
    tcode = np.mintypecode([a.dtype.char for a in arrays])
    out = _valarray(np.shape(arrays[0]), value=fillvalue, typecode=tcode)
    np.place(out, cond, f(*temp))
    if f2 is not None:
        temp = tuple(np.extract(~cond, arr) for arr in arrays)
        np.place(out, ~cond, f2(*temp))

    return out


# Work around for complex chnges in gammaln in 1.0.0.  loggamma introduced in 0.18. 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def compute_KLD(self, tup, batch_sz):
        kappa = tup['kappa']
        d = self.lat_dim

        rt_bag = []
        # const = torch.log(torch.tensor(3.1415926)) * d / 2 + torch.log(torch.tensor(2.0)) \
        #         - torch.tensor(sp.loggamma(d / 2).real) - (d / 2) * torch.log(torch.tensor(2 * 3.1415926))

        const = torch.tensor(
            np.log(np.pi) * d / 2 + np.log(2) - sp.loggamma(d / 2).real - (d / 2) * np.log(2 * np.pi)).to(
            device)
        d = torch.tensor([d], dtype=torch.float).to(device)
        batchsz = kappa.size()[0]

        rt_tensor = torch.zeros(batchsz)
        for k_idx in range(batchsz):
            k = kappa[k_idx]
            # print(k)
            # print(k)
            # print(d)
            first = k * bessel_iv(d / 2, k) / bessel_iv(d / 2 - 1, k)
            second = (d / 2 - 1) * torch.log(k) - torch.log(bessel_iv(d / 2 - 1, k))
            combin = first + second + const
            rt_tensor[k_idx] = combin
            # rt_bag.append(combin)
        return rt_tensor.to(device)
        # return torch.tensor(rt_bag,requires_grad=True).to(device) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def KL_guu(k, d):
    kld = k * ((sp.iv(d / 2.0 + 1.0, k) \
                + sp.iv(d / 2.0, k) * d / (2.0 * k)) / sp.iv(d / 2.0, k) - d / (2.0 * k)) \
          + d * np.log(k) / 2.0 - np.log(sp.iv(d / 2.0, k)) \
          - sp.loggamma(d / 2 + 1) - d * np.log(2) / 2

    return kld 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def _vmfKL(k, d):
    return k * ((sp.iv(d / 2.0 + 1.0, k) \
                 + sp.iv(d / 2.0, k) * d / (2.0 * k)) / sp.iv(d / 2.0, k) - d / (2.0 * k)) \
           + d * np.log(k) / 2.0 - np.log(sp.iv(d / 2.0, k)) \
           - sp.loggamma(d / 2 + 1) - d * np.log(2) / 2 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def _vmf_kld(k, d):
        tmp = (k * ((sp.iv(d / 2.0 + 1.0, k) + sp.iv(d / 2.0, k) * d / (2.0 * k)) / sp.iv(d / 2.0, k) - d / (2.0 * k)) \
               + d * np.log(k) / 2.0 - np.log(sp.iv(d / 2.0, k)) \
               - sp.loggamma(d / 2 + 1) - d * np.log(2) / 2).real
        if tmp != tmp:
            exit()
        return np.array([tmp]) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def _vmf_kld_davidson(k, d):
        """
        This should be the correct KLD.
        Empirically we find that _vmf_kld (as in the Guu paper) only deviates a little (<2%) in most cases we use.
        """
        tmp = k * sp.iv(d / 2, k) / sp.iv(d / 2 - 1, k) + (d / 2 - 1) * torch.log(k) - torch.log(
            sp.iv(d / 2 - 1, k)) + np.log(np.pi) * d / 2 + np.log(2) - sp.loggamma(d / 2).real - (d / 2) * np.log(
            2 * np.pi)
        if tmp != tmp:
            exit()
        return np.array([tmp]) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def _vmf_kld(k, d):
        tmp = (k * ((sp.iv(d / 2.0 + 1.0, k) + sp.iv(d / 2.0, k) * d / (2.0 * k)) / sp.iv(d / 2.0, k) - d / (2.0 * k)) \
               + d * np.log(k) / 2.0 - np.log(sp.iv(d / 2.0, k)) \
               - sp.loggamma(d / 2 + 1) - d * np.log(2) / 2).real
        return tmp 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def test_identities1():
    # test the identity exp(loggamma(z)) = gamma(z)
    x = np.array([-99.5, -9.5, -0.5, 0.5, 9.5, 99.5])
    y = x.copy()
    x, y = np.meshgrid(x, y)
    z = (x + 1J*y).flatten()
    dataset = np.vstack((z, gamma(z))).T

    def f(z):
        return np.exp(loggamma(z))

    FuncData(f, dataset, 0, 1, rtol=1e-14, atol=1e-14).check() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def test_identities2():
    # test the identity loggamma(z + 1) = log(z) + loggamma(z)
    x = np.array([-99.5, -9.5, -0.5, 0.5, 9.5, 99.5])
    y = x.copy()
    x, y = np.meshgrid(x, y)
    z = (x + 1J*y).flatten()
    dataset = np.vstack((z, np.log(z) + loggamma(z))).T

    def f(z):
        return loggamma(z + 1)

    FuncData(f, dataset, 0, 1, rtol=1e-14, atol=1e-14).check() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def test_complex_dispatch_realpart():
    # Test that the real parts of loggamma and gammaln agree on the
    # real axis.
    x = np.r_[-np.logspace(10, -10), np.logspace(-10, 10)] + 0.5

    dataset = np.vstack((x, gammaln(x))).T

    def f(z):
        z = np.array(z, dtype='complex128')
        return loggamma(z).real

    FuncData(f, dataset, 0, 1, rtol=1e-14, atol=1e-14).check() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import loggamma [as 别名]
def test_real_dispatch():
    x = np.logspace(-10, 10) + 0.5
    dataset = np.vstack((x, gammaln(x))).T

    FuncData(loggamma, dataset, 0, 1, rtol=1e-14, atol=1e-14).check()
    assert_(loggamma(0) == np.inf)
    assert_(np.isnan(loggamma(-1)))