# Python scipy.special.digamma() Examples

The following are 30 code examples of scipy.special.digamma(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module , or try the search function .
Example #1
```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) ```
Example #2
```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) ```
Example #3
```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
# ------------------------------------------------------------------------------ ```
Example #4
```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() ```
Example #5
```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() ```
Example #6
```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
# ------------------------------------------------------------------------------ ```
Example #7
```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() ```
Example #8
 Source File: entropy_estimators.py    From scikit-feature with GNU General Public License v2.0 6 votes
```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) ```
Example #9
 Source File: entropy_estimators.py    From scikit-feature with GNU General Public License v2.0 6 votes
```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) ```
Example #10
 Source File: negative_binomial.py    From DiscoSnp with GNU Affero General Public License v3.0 6 votes
```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))

Example #11
```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) ```
Example #12
```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) ```
Example #13
```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
######################################################### ```
Example #14
```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) ```
Example #15
```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 ```
Example #16
```def _digammainv(y):
# Inverse of the digamma function (real positive arguments only).
# This function is used in the `fit` method of `gamma_gen`.
# The function uses either optimize.fsolve or optimize.newton
# to solve `sc.digamma(x) - y = 0`.  There is probably room for
# improvement, but currently it works over a wide range of y:
#    >>> y = 64*np.random.randn(1000000)
#    >>> y.min(), y.max()
#    (-311.43592651416662, 351.77388222276869)
#    x = [_digammainv(t) for t in y]
#    np.abs(sc.digamma(x) - y).max()
#    1.1368683772161603e-13
#
_em = 0.5772156649015328606065120
func = lambda x: sc.digamma(x) - y
if y > -0.125:
x0 = np.exp(y) + 0.5
if y < 10:
# Some experimentation shows that newton reliably converges
# must faster than fsolve in this y range.  For larger y,
# newton sometimes fails to converge.
value = optimize.newton(func, x0, tol=1e-10)
return value
elif y > -3:
x0 = np.exp(y/2.332) + 0.08661
else:
x0 = 1.0 / (-y - _em)

value, info, ier, mesg = optimize.fsolve(func, x0, xtol=1e-11,
full_output=True)
if ier != 1:
raise RuntimeError("_digammainv: fsolve failed, y = %r" % y)

return value[0]

## Gamma (Use MATLAB and MATHEMATICA (b=theta=scale, a=alpha=shape) definition)

## gamma(a, loc, scale)  with a an integer is the Erlang distribution
## gamma(1, loc, scale)  is the Exponential distribution
## gamma(df/2, 0, 2) is the chi2 distribution with df degrees of freedom. ```
Example #17
```def _stats(self, c):
# See, for example, "A Statistical Study of Log-Gamma Distribution", by
# Ping Shing Chan (thesis, McMaster University, 1993).
mean = sc.digamma(c)
var = sc.polygamma(1, c)
skewness = sc.polygamma(2, c) / np.power(var, 1.5)
excess_kurtosis = sc.polygamma(3, c) / (var*var)
return mean, var, skewness, excess_kurtosis ```
Example #18
```def _score_nbin(self, params, Q=0):
"""
Score vector for NB2 model
"""
if self._transparams: # lnalpha came in during fit
alpha = np.exp(params[-1])
else:
alpha = params[-1]
params = params[:-1]
exog = self.exog
y = self.endog[:,None]
mu = self.predict(params)[:,None]
a1 = 1/alpha * mu**Q
prob = a1 / (a1 + mu)  # a1 aka "size" in _ll_nbin
if Q == 1:  # nb1
# Q == 1 --> a1 = mu / alpha --> prob = 1 / (alpha + 1)
dgpart = digamma(y + a1) - digamma(a1)
dparams = exog * a1 * (np.log(prob) +
dgpart)
dalpha = ((alpha * (y - mu * np.log(prob) -
mu*(dgpart + 1)) -
mu * (np.log(prob) +
dgpart))/
(alpha**2*(alpha + 1))).sum()

elif Q == 0:  # nb2
dgpart = digamma(y + a1) - digamma(a1)
dparams = exog*a1 * (y-mu)/(mu+a1)
da1 = -alpha**-2
dalpha = (dgpart + np.log(a1)
- np.log(a1+mu) - (y-mu)/(a1+mu)).sum() * da1

#multiply above by constant outside sum to reduce rounding error
if self._transparams:
return np.r_[dparams.sum(0), dalpha*alpha]
else:
return np.r_[dparams.sum(0), dalpha] ```
Example #19
```def test_fix_fit_gamma(self):
x = np.arange(1, 6)
meanlog = np.log(x).mean()

# A basic test of gamma.fit with floc=0.
floc = 0
a, loc, scale = stats.gamma.fit(x, floc=floc)
s = np.log(x.mean()) - meanlog
assert_almost_equal(np.log(a) - special.digamma(a), s, decimal=5)
assert_equal(loc, floc)
assert_almost_equal(scale, x.mean()/a, decimal=8)

# Regression tests for gh-2514.
# The problem was that if `floc=0` was given, any other fixed
# parameters were ignored.
f0 = 1
floc = 0
a, loc, scale = stats.gamma.fit(x, f0=f0, floc=floc)
assert_equal(a, f0)
assert_equal(loc, floc)
assert_almost_equal(scale, x.mean()/a, decimal=8)

f0 = 2
floc = 0
a, loc, scale = stats.gamma.fit(x, f0=f0, floc=floc)
assert_equal(a, f0)
assert_equal(loc, floc)
assert_almost_equal(scale, x.mean()/a, decimal=8)

# loc and scale fixed.
floc = 0
fscale = 2
a, loc, scale = stats.gamma.fit(x, floc=floc, fscale=fscale)
assert_equal(loc, floc)
assert_equal(scale, fscale)
c = meanlog - np.log(fscale)
assert_almost_equal(special.digamma(a), c) ```
Example #20
```def test_digamma(self):
assert_mpmath_equal(sc.digamma,
_exception_to_nan(mpmath.digamma),
[Arg()],
dps=50) ```
Example #21
```def test_digamma_complex(self):
assert_mpmath_equal(sc.digamma,
_time_limited()(_exception_to_nan(mpmath.digamma)),
[ComplexArg()],
n=200) ```
Example #22
```def _estimate_log_weights(self):
if self.weight_concentration_prior_type == 'dirichlet_process':
digamma_sum = digamma(self.weight_concentration_[0] +
self.weight_concentration_[1])
digamma_a = digamma(self.weight_concentration_[0])
digamma_b = digamma(self.weight_concentration_[1])
return (digamma_a - digamma_sum +
np.hstack((0, np.cumsum(digamma_b - digamma_sum)[:-1])))
else:
# case Variationnal Gaussian mixture with dirichlet distribution
return (digamma(self.weight_concentration_) -
digamma(np.sum(self.weight_concentration_))) ```
Example #23
```def _estimate_log_prob(self, X):
_, n_features = X.shape
# We remove `n_features * np.log(self.degrees_of_freedom_)` because
# the precision matrix is normalized
log_gauss = (_estimate_log_gaussian_prob(
X, self.means_, self.precisions_cholesky_, self.covariance_type) -
.5 * n_features * np.log(self.degrees_of_freedom_))

log_lambda = n_features * np.log(2.) + np.sum(digamma(
.5 * (self.degrees_of_freedom_ -
np.arange(0, n_features)[:, np.newaxis])), 0)

return log_gauss + .5 * (log_lambda -
n_features / self.mean_precision_) ```
Example #24
```def _digamma_cpu(x, dtype):
from scipy import special
return numpy.vectorize(special.digamma, otypes=[dtype])(x) ```
Example #25
```def label(self):
return 'digamma' ```
Example #26
```def forward_cpu(self, x):
global _digamma_cpu
if _digamma_cpu is None:
try:
from scipy import special
_digamma_cpu = special.digamma
except ImportError:
raise ImportError('SciPy is not available. Forward computation'
' of digamma can not be done.')
self.retain_inputs((0,))
return utils.force_array(_digamma_cpu(x[0]), dtype=x[0].dtype), ```
Example #27
```def forward_gpu(self, x):
self.retain_inputs((0,))
return utils.force_array(
cuda.cupyx.scipy.special.digamma(x[0]), dtype=x[0].dtype), ```
Example #28
```def digamma(x):
"""Digamma function.

.. note::
Forward computation in CPU can not be done if
`SciPy <https://www.scipy.org/>`_ is not available.

Args:
x (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.

Returns:
~chainer.Variable: Output variable.
"""
return DiGamma().apply((x,))[0] ```
Example #29
 Source File: lda.py    From numpy-ml with GNU General Public License v3.0 5 votes
```def _maximize_alpha(self, max_iters=1000, tol=0.1):
"""
Optimize alpha using Blei's O(n) Newton-Raphson modification
for a Hessian with special structure
"""
D = self.D
T = self.T

alpha = self.alpha
gamma = self.gamma

for _ in range(max_iters):
alpha_old = alpha

g = D * (digamma(np.sum(alpha)) - digamma(alpha)) + np.sum(
digamma(gamma) - np.tile(digamma(np.sum(gamma, axis=1)), (T, 1)).T,
axis=0,
)

#  Calculate Hessian diagonal component
h = -D * polygamma(1, alpha)

#  Calculate Hessian constant component
z = D * polygamma(1, np.sum(alpha))

#  Calculate constant
c = np.sum(g / h) / (z ** (-1.0) + np.sum(h ** (-1.0)))

#  Update alpha
alpha = alpha - (g - c) / h

#  Check convergence
if np.sqrt(np.mean(np.square(alpha - alpha_old))) < tol:
break

return alpha ```
Example #30
 Source File: lda.py    From numpy-ml with GNU General Public License v3.0 5 votes
```def dg(gamma, d, t):
"""
E[log X_t] where X_t ~ Dir
"""
return digamma(gamma[d, t]) - digamma(np.sum(gamma[d, :])) ```