Python special.digamma方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def mi(x, y, k=3, base=2):
  """Mutual information of x and y.
  x,y should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
  if x is a one-dimensional scalar and we have four samples.
  """
  assert len(x)==len(y), 'Lists should have same length.'
  assert k <= len(x) - 1, 'Set k smaller than num samples - 1.'
  intens = 1e-10 # Small noise to break degeneracy, see doc.
  x = [list(p + intens*nr.rand(len(x[0]))) for p in x]
  y = [list(p + intens*nr.rand(len(y[0]))) for p in y]
  points = zip2(x,y)
  # Find nearest neighbors in joint space, p=inf means max-norm.
  tree = ss.cKDTree(points)
  dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points]
  a = avgdigamma(x,dvec)
  b = avgdigamma(y,dvec)
  c = digamma(k)
  d = digamma(len(x))
  return (-a-b+c+d) / log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def cmi(x, y, z, k=3, base=2):
  """Mutual information of x and y, conditioned on z
  x,y,z should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
  if x is a one-dimensional scalar and we have four samples
  """
  assert len(x)==len(y), 'Lists should have same length.'
  assert k <= len(x) - 1, 'Set k smaller than num samples - 1.'
  intens = 1e-10 # Small noise to break degeneracy, see doc.
  x = [list(p + intens*nr.rand(len(x[0]))) for p in x]
  y = [list(p + intens*nr.rand(len(y[0]))) for p in y]
  z = [list(p + intens*nr.rand(len(z[0]))) for p in z]
  points = zip2(x,y,z)
  # Find nearest neighbors in joint space, p=inf means max-norm.
  tree = ss.cKDTree(points)
  dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points]
  a = avgdigamma(zip2(x,z), dvec)
  b = avgdigamma(zip2(y,z), dvec)
  c = avgdigamma(z,dvec)
  d = digamma(k)
  return (-a-b+c+d) / log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def test_rgamma_zeros():
    # Test around the zeros at z = 0, -1, -2, ...,  -169. (After -169 we
    # get values that are out of floating point range even when we're
    # within 0.1 of the zero.)

    # Can't use too many points here or the test takes forever.
    dx = np.r_[-np.logspace(-1, -13, 3), 0, np.logspace(-13, -1, 3)]
    dy = dx.copy()
    dx, dy = np.meshgrid(dx, dy)
    dz = dx + 1j*dy
    zeros = np.arange(0, -170, -1).reshape(1, 1, -1)
    z = (zeros + np.dstack((dz,)*zeros.size)).flatten()
    dataset = []
    with mpmath.workdps(100):
        for z0 in z:
            dataset.append((z0, complex(mpmath.rgamma(z0))))

    dataset = np.array(dataset)
    FuncData(sc.rgamma, dataset, 0, 1, rtol=1e-12).check()


# ------------------------------------------------------------------------------
# digamma
# ------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def test_digamma_roots():
    # Test the special-cased roots for digamma.
    root = mpmath.findroot(mpmath.digamma, 1.5)
    roots = [float(root)]
    root = mpmath.findroot(mpmath.digamma, -0.5)
    roots.append(float(root))
    roots = np.array(roots)

    # If we test beyond a radius of 0.24 mpmath will take forever.
    dx = np.r_[-0.24, -np.logspace(-1, -15, 10), 0, np.logspace(-15, -1, 10), 0.24]
    dy = dx.copy()
    dx, dy = np.meshgrid(dx, dy)
    dz = dx + 1j*dy
    z = (roots + np.dstack((dz,)*roots.size)).flatten()
    dataset = []
    with mpmath.workdps(30):
        for z0 in z:
            dataset.append((z0, complex(mpmath.digamma(z0))))

    dataset = np.array(dataset)
    FuncData(sc.digamma, dataset, 0, 1, rtol=1e-14).check() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def test_digamma_negreal():
    # Test digamma around the negative real axis. Don't do this in
    # TestSystematic because the points need some jiggering so that
    # mpmath doesn't take forever.

    digamma = exception_to_nan(mpmath.digamma)

    x = -np.logspace(300, -30, 100)
    y = np.r_[-np.logspace(0, -3, 5), 0, np.logspace(-3, 0, 5)]
    x, y = np.meshgrid(x, y)
    z = (x + 1j*y).flatten()

    dataset = []
    with mpmath.workdps(40):
        for z0 in z:
            res = digamma(z0)
            dataset.append((z0, complex(res)))
    dataset = np.asarray(dataset)

    FuncData(sc.digamma, dataset, 0, 1, rtol=1e-13).check() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def test_digamma_boundary():
    # Check that there isn't a jump in accuracy when we switch from
    # using the asymptotic series to the reflection formula.

    x = -np.logspace(300, -30, 100)
    y = np.array([-6.1, -5.9, 5.9, 6.1])
    x, y = np.meshgrid(x, y)
    z = (x + 1j*y).flatten()

    dataset = []
    with mpmath.workdps(30):
        for z0 in z:
            res = mpmath.digamma(z0)
            dataset.append((z0, complex(res)))
    dataset = np.asarray(dataset)

    FuncData(sc.digamma, dataset, 0, 1, rtol=1e-13).check()


# ------------------------------------------------------------------------------
# gammainc
# ------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def avgdigamma(data, dvec, leaf_size=16):
    """Convenience function for finding expectation value of <psi(nx)> given
    some number of neighbors in some radius in a marginal space.

    Parameters
    ----------
    points : numpy.ndarray
    dvec : array_like (n_points,)
    Returns
    -------
    avgdigamma : float
        expectation value of <psi(nx)>
    """
    tree = BallTree(data, leaf_size=leaf_size, p=float('inf'))

    n_points = tree.query_radius(data, dvec - EPS, count_only=True)

    return digamma(n_points).mean() 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def mi(x, y, k=3, base=2):
    """
    Mutual information of x and y; x, y should be a list of vectors, e.g. x = [[1.3],[3.7],[5.1],[2.4]]
    if x is a one-dimensional scalar and we have four samples
    """

    assert len(x) == len(y), "Lists should have same length"
    assert k <= len(x) - 1, "Set k smaller than num. samples - 1"
    intens = 1e-10  # small noise to break degeneracy, see doc.
    x = [list(p + intens * nr.rand(len(x[0]))) for p in x]
    y = [list(p + intens * nr.rand(len(y[0]))) for p in y]
    points = zip2(x, y)
    # Find nearest neighbors in joint space, p=inf means max-norm
    tree = ss.cKDTree(points)
    dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points]
    a, b, c, d = avgdigamma(x, dvec), avgdigamma(y, dvec), digamma(k), digamma(len(x))
    return (-a-b+c+d)/log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def cmi(x, y, z, k=3, base=2):
    """
    Mutual information of x and y, conditioned on z; x, y, z should be a list of vectors, e.g. x = [[1.3],[3.7],[5.1],[2.4]]
    if x is a one-dimensional scalar and we have four samples
    """

    assert len(x) == len(y), "Lists should have same length"
    assert k <= len(x) - 1, "Set k smaller than num. samples - 1"
    intens = 1e-10  # small noise to break degeneracy, see doc.
    x = [list(p + intens * nr.rand(len(x[0]))) for p in x]
    y = [list(p + intens * nr.rand(len(y[0]))) for p in y]
    z = [list(p + intens * nr.rand(len(z[0]))) for p in z]
    points = zip2(x, y, z)
    # Find nearest neighbors in joint space, p=inf means max-norm
    tree = ss.cKDTree(points)
    dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points]
    a, b, c, d = avgdigamma(zip2(x, z), dvec), avgdigamma(zip2(y, z), dvec), avgdigamma(z, dvec), digamma(k)
    return (-a-b+c+d)/log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def r_derv(r_var, vec):
    ''' Function that represents the derivative of the neg bin likelihood wrt r
    @param r: The value of r in the derivative of the likelihood wrt r
    @param vec: The data vector used in the likelihood
    '''
    if not r_var or not vec:
        raise ValueError("r parameter and data must be specified")

    if r_var <= 0:
        raise ValueError("r must be strictly greater than 0")

    total_sum = 0
    obs_mean = np.mean(vec)  # Save the mean of the data
    n_pop = float(len(vec))  # Save the length of the vector, n_pop

    for obs in vec:
        total_sum += digamma(obs + r_var)

    total_sum -= n_pop*digamma(r_var)
    total_sum += n_pop*math.log(r_var / (r_var + obs_mean))

    return total_sum 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def mutual_info(x, y, k=3, base=2):
        """ Mutual information of x and y
            x, y should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
            if x is a one-dimensional scalar and we have four samples
        """
        x=entropyc.__reshape(x)
        y=entropyc.__reshape(y)
        assert len(x) == len(y), "Lists should have same length"
        assert k <= len(x) - 1, "Set k smaller than num. samples - 1"
        intens = 1e-10  # small noise to break degeneracy, see doc.
        x = [list(p + intens * nr.rand(len(x[0]))) for p in x]
        y = [list(p + intens * nr.rand(len(y[0]))) for p in y]
        points = zip2(x, y)
        # Find nearest neighbors in joint space, p=inf means max-norm
        tree = ss.cKDTree(points)
        dvec = [tree.query(point, k + 1, p=float('inf'))[0][k] for point in points]
        a, b, c, d = avgdigamma(x, dvec), avgdigamma(y, dvec), digamma(k), digamma(len(x))
        return (-a - b + c + d) / log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def cond_mutual_info(x, y, z, k=3, base=2):
        """ Mutual information of x and y, conditioned on z
            x, y, z should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
            if x is a one-dimensional scalar and we have four samples
        """
        x=entropyc.__reshape(x)
        y=entropyc.__reshape(y)
        z=entropyc.__reshape(z)
        assert len(x) == len(y), "Lists should have same length"
        assert k <= len(x) - 1, "Set k smaller than num. samples - 1"
        intens = 1e-10  # small noise to break degeneracy, see doc.
        x = [list(p + intens * nr.rand(len(x[0]))) for p in x]
        y = [list(p + intens * nr.rand(len(y[0]))) for p in y]
        z = [list(p + intens * nr.rand(len(z[0]))) for p in z]
        points = zip2(x, y, z)
        # Find nearest neighbors in joint space, p=inf means max-norm
        tree = ss.cKDTree(points)
        dvec = [tree.query(point, k + 1, p=float('inf'))[0][k] for point in points]
        a, b, c, d = avgdigamma(zip2(x, z), dvec), avgdigamma(zip2(y, z), dvec), avgdigamma(z, dvec), digamma(k)
        return (-a - b + c + d) / log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def calc_mergeLP(self, LP, kA, kB):
    ''' Calculate and return the new local parameters for a merged configuration
          that combines topics kA and kB

        Returns
        ---------
        LP : dict of local params, with updated fields and only K-1 comps
    '''
    K = LP['DocTopicCount'].shape[1]
    assert kA < kB
    LP = copy.deepcopy(LP)
    for key in ['DocTopicCount', 'word_variational', 'alphaPi']:
      LP[key][:,kA] = LP[key][:,kA] + LP[key][:,kB]
      LP[key] = np.delete(LP[key], kB, axis=1)
    LP['E_logPi'] = digamma(LP['alphaPi']) \
                    - digamma(LP['alphaPi'].sum(axis=1))[:,np.newaxis]
    assert LP['word_variational'].shape[1] == K-1
    return LP


  ######################################################### Verify Z terms
  ######################################################### 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def entropy(x, k=3, base=2):
  """The classic K-L k-nearest neighbor continuous entropy estimator.
  x should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
  if x is a one-dimensional scalar and we have four samples
  """
  assert k <= len(x)-1, 'Set k smaller than num. samples - 1.'
  d = len(x[0])
  N = len(x)
  intens = 1e-10 # Small noise to break degeneracy, see doc.
  x = [list(p + intens*nr.rand(len(x[0]))) for p in x]
  tree = ss.cKDTree(x)
  nn = [tree.query(point, k+1, p=float('inf'))[0][k] for point in x]
  const = digamma(N) - digamma(k) + d*log(2)
  return (const + d*np.mean(map(log, nn))) / log(base) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import digamma [as 别名]
def avgdigamma(points, dvec):
  # This part finds number of neighbors in some radius in the marginal space
  # returns expectation value of <psi(nx)>.
  N = len(points)
  tree = ss.cKDTree(points)
  avg = 0.
  for i in xrange(N):
    dist = dvec[i]
    # Subtlety, we don't include the boundary point,
    # but we are implicitly adding 1 to kraskov def bc center point is included.
    num_points = len(tree.query_ball_point(points[i], dist-1e-15,
      p=float('inf')))
    avg += digamma(num_points) / N
  return avg