pyerrors.special
1import scipy 2import numpy as np 3from autograd.extend import primitive, defvjp 4from autograd.scipy.special import j0, y0, j1, y1, jn, yn, i0, i1, iv, ive, beta, betainc, betaln 5from autograd.scipy.special import polygamma, psi, digamma, gamma, gammaln, gammainc, gammaincc, gammasgn, rgamma, multigammaln 6from autograd.scipy.special import erf, erfc, erfinv, erfcinv, logit, expit, logsumexp 7 8 9__all__ = ["beta", "betainc", "betaln", 10 "polygamma", "psi", "digamma", "gamma", "gammaln", "gammainc", "gammaincc", "gammasgn", "rgamma", "multigammaln", 11 "kn", "j0", "y0", "j1", "y1", "jn", "yn", "i0", "i1", "iv", "ive", 12 "erf", "erfc", "erfinv", "erfcinv", "logit", "expit", "logsumexp"] 13 14 15@primitive 16def kn(n, x): 17 """Modified Bessel function of the second kind of integer order n""" 18 if int(n) != n: 19 raise TypeError("The order 'n' needs to be an integer.") 20 return scipy.special.kn(n, x) 21 22 23defvjp(kn, None, lambda ans, n, x: lambda g: - g * 0.5 * (kn(np.abs(n - 1), x) + kn(n + 1, x)))
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
beta(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
beta(a, b, out=None)
Beta function.
This function is defined in 1 as
$$B(a, b) = \int_0^1 t^{a-1}(1-t)^{b-1}dt = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)},$$
where \( \Gamma \) is the gamma function.
Parameters
- a, b (array_like): Real-valued arguments
- out (ndarray, optional): Optional output array for the function result
Returns
- scalar or ndarray: Value of the beta function
See Also
gamma
: the gamma function
betainc
: the regularized incomplete beta function
betaln
: the natural logarithm of the absolute
value of the beta function
References
Examples
>>> import scipy.special as sc
The beta function relates to the gamma function by the definition given above:
>>> sc.beta(2, 3)
0.08333333333333333
>>> sc.gamma(2)*sc.gamma(3)/sc.gamma(2 + 3)
0.08333333333333333
As this relationship demonstrates, the beta function is symmetric:
>>> sc.beta(1.7, 2.4)
0.16567527689031739
>>> sc.beta(2.4, 1.7)
0.16567527689031739
This function satisfies \( B(1, b) = 1/b \):
>>> sc.beta(1, 4)
0.25
-
NIST Digital Library of Mathematical Functions, Eq. 5.12.1. https://dlmf.nist.gov/5.12 ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
betainc(x1, x2, x3, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
betainc(a, b, x, out=None)
Regularized incomplete beta function.
Computes the regularized incomplete beta function, defined as 1:
$$I_x(a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int_0^x t^{a-1}(1-t)^{b-1}dt,$$
for \( 0 \leq x \leq 1 \).
This function is the cumulative distribution function for the beta distribution; its range is [0, 1].
Parameters
- a, b (array_like): Positive, real-valued parameters
- x (array_like): Real-valued such that \( 0 \leq x \leq 1 \), the upper limit of integration
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Value of the regularized incomplete beta function
See Also
beta
: beta function
betaincinv
: inverse of the regularized incomplete beta function
betaincc
: complement of the regularized incomplete beta function
scipy.stats.beta
: beta distribution
Notes
The term regularized in the name of this function refers to the
scaling of the function by the gamma function terms shown in the
formula. When not qualified as regularized, the name incomplete
beta function often refers to just the integral expression,
without the gamma terms. One can use the function beta
from
scipy.special
to get this "nonregularized" incomplete beta
function by multiplying the result of betainc(a, b, x)
by
beta(a, b)
.
References
Examples
Let \( B(a, b) \) be the beta
function.
>>> import scipy.special as sc
The coefficient in terms of gamma
is equal to
\( 1/B(a, b) \). Also, when \( x=1 \)
the integral is equal to \( B(a, b) \).
Therefore, \( I_{x=1}(a, b) = 1 \) for any \( a, b \).
>>> sc.betainc(0.2, 3.5, 1.0)
1.0
It satisfies
\( I_x(a, b) = x^a F(a, 1-b, a+1, x)/ (aB(a, b)) \),
where \( F \) is the hypergeometric function hyp2f1
:
>>> a, b, x = 1.4, 3.1, 0.5
>>> x**a * sc.hyp2f1(a, 1 - b, a + 1, x)/(a * sc.beta(a, b))
0.8148904036225295
>>> sc.betainc(a, b, x)
0.8148904036225296
This functions satisfies the relationship \( I_x(a, b) = 1 - I_{1-x}(b, a) \):
>>> sc.betainc(2.2, 3.1, 0.4)
0.49339638807619446
>>> 1 - sc.betainc(3.1, 2.2, 1 - 0.4)
0.49339638807619446
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/8.17 ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
betaln(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
betaln(a, b, out=None)
Natural logarithm of absolute value of beta function.
Computes ln(abs(beta(a, b)))
.
Parameters
- a, b (array_like): Positive, real-valued parameters
- out (ndarray, optional): Optional output array for function values
Returns
- scalar or ndarray: Value of the betaln function
See Also
gamma
: the gamma function
betainc
: the regularized incomplete beta function
beta
: the beta function
Examples
>>> import numpy as np
>>> from scipy.special import betaln, beta
Verify that, for moderate values of a
and b
, betaln(a, b)
is the same as log(beta(a, b))
:
>>> betaln(3, 4)
-4.0943445622221
>>> np.log(beta(3, 4))
-4.0943445622221
In the following beta(a, b)
underflows to 0, so we can't compute
the logarithm of the actual value.
>>> a = 400
>>> b = 900
>>> beta(a, b)
0.0
We can compute the logarithm of beta(a, b)
by using betaln
:
>>> betaln(a, b)
-804.3069951764146
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
Polygamma functions.
Defined as \( \psi^{(n)}(x) \) where \( \psi \) is the
digamma
function. See [dlmf]_ for details.
Parameters
- n (array_like): The order of the derivative of the digamma function; must be integral
- x (array_like): Real valued input
Returns
- ndarray: Function results
See Also
References
.. [dlmf] NIST, Digital Library of Mathematical Functions, https://dlmf.nist.gov/5.15
Examples
>>> from scipy import special
>>> x = [2, 3, 25.5]
>>> special.polygamma(1, x)
array([ 0.64493407, 0.39493407, 0.03999467])
>>> special.polygamma(0, x) == special.psi(x)
array([ True, True, True], dtype=bool)
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
psi(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
psi(z, out=None)
The digamma function.
The logarithmic derivative of the gamma function evaluated at z
.
Parameters
- z (array_like): Real or complex argument.
- out (ndarray, optional):
Array for the computed values of
psi
.
Returns
- digamma (scalar or ndarray):
Computed values of
psi
.
Notes
For large values not close to the negative real axis, psi
is
computed using the asymptotic series (5.11.2) from 1. For small
arguments not close to the negative real axis, the recurrence
relation (5.5.2) from 2 is used until the argument is large
enough to use the asymptotic series. For values close to the
negative real axis, the reflection formula (5.5.4) from 3 is
used first. Note that psi
has a family of zeros on the
negative real axis which occur between the poles at nonpositive
integers. Around the zeros the reflection formula suffers from
cancellation and the implementation loses precision. The sole
positive zero and the first negative zero, however, are handled
separately by precomputing series expansions using 4, so the
function should maintain full accuracy around the origin.
References
Examples
>>> from scipy.special import psi
>>> z = 3 + 4j
>>> psi(z)
(1.55035981733341+1.0105022091860445j)
Verify psi(z) = psi(z + 1) - 1/z:
>>> psi(z + 1) - 1/z
(1.55035981733341+1.0105022091860445j)
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩
-
Fredrik Johansson and others. "mpmath: a Python library for arbitrary-precision floating-point arithmetic" (Version 0.19) http://mpmath.org/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
psi(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
psi(z, out=None)
The digamma function.
The logarithmic derivative of the gamma function evaluated at z
.
Parameters
- z (array_like): Real or complex argument.
- out (ndarray, optional):
Array for the computed values of
psi
.
Returns
- digamma (scalar or ndarray):
Computed values of
psi
.
Notes
For large values not close to the negative real axis, psi
is
computed using the asymptotic series (5.11.2) from 1. For small
arguments not close to the negative real axis, the recurrence
relation (5.5.2) from 2 is used until the argument is large
enough to use the asymptotic series. For values close to the
negative real axis, the reflection formula (5.5.4) from 3 is
used first. Note that psi
has a family of zeros on the
negative real axis which occur between the poles at nonpositive
integers. Around the zeros the reflection formula suffers from
cancellation and the implementation loses precision. The sole
positive zero and the first negative zero, however, are handled
separately by precomputing series expansions using 4, so the
function should maintain full accuracy around the origin.
References
Examples
>>> from scipy.special import psi
>>> z = 3 + 4j
>>> psi(z)
(1.55035981733341+1.0105022091860445j)
Verify psi(z) = psi(z + 1) - 1/z:
>>> psi(z + 1) - 1/z
(1.55035981733341+1.0105022091860445j)
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩
-
NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩
-
Fredrik Johansson and others. "mpmath: a Python library for arbitrary-precision floating-point arithmetic" (Version 0.19) http://mpmath.org/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
gamma(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
gamma(z, out=None)
gamma function.
The gamma function is defined as
$$\Gamma(z) = \int_0^\infty t^{z-1} e^{-t} dt$$
for \( \Re(z) > 0 \) and is extended to the rest of the complex plane by analytic continuation. See [dlmf]_ for more details.
Parameters
- z (array_like): Real or complex valued argument
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Values of the gamma function
Notes
The gamma function is often referred to as the generalized factorial since \( \Gamma(n + 1) = n! \) for natural numbers \( n \). More generally it satisfies the recurrence relation \( \Gamma(z + 1) = z \cdot \Gamma(z) \) for complex \( z \), which, combined with the fact that \( \Gamma(1) = 1 \), implies the above identity for \( z = n \).
References
.. [dlmf] NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5.2#E1
Examples
>>> import numpy as np
>>> from scipy.special import gamma, factorial
>>> gamma([0, 0.5, 1, 5])
array([ inf, 1.77245385, 1. , 24. ])
>>> z = 2.5 + 1j
>>> gamma(z)
(0.77476210455108352+0.70763120437959293j)
>>> gamma(z+1), z*gamma(z) # Recurrence property
((1.2292740569981171+2.5438401155000685j),
(1.2292740569981158+2.5438401155000658j))
>>> gamma(0.5)**2 # gamma(0.5) = sqrt(pi)
3.1415926535897927
Plot gamma(x) for real x
>>> x = np.linspace(-3.5, 5.5, 2251)
>>> y = gamma(x)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'b', alpha=0.6, label='gamma(x)')
>>> k = np.arange(1, 7)
>>> plt.plot(k, factorial(k-1), 'k*', alpha=0.6,
... label='(x-1)!, x = 1, 2, ...')
>>> plt.xlim(-3.5, 5.5)
>>> plt.ylim(-10, 25)
>>> plt.grid()
>>> plt.xlabel('x')
>>> plt.legend(loc='lower right')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
gammaln(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
gammaln(x, out=None)
Logarithm of the absolute value of the gamma function.
Defined as
$$\ln(\lvert\Gamma(x)\rvert)$$
where \( \Gamma \) is the gamma function. For more details on the gamma function, see [dlmf]_.
Parameters
- x (array_like): Real argument
- out (ndarray, optional): Optional output array for the function results
Returns
- scalar or ndarray: Values of the log of the absolute value of gamma
See Also
gammasgn
: sign of the gamma function
loggamma
: principal branch of the logarithm of the gamma function
Notes
It is the same function as the Python standard library function
math.lgamma()
.
When used in conjunction with gammasgn
, this function is useful
for working in logspace on the real axis without having to deal
with complex numbers via the relation exp(gammaln(x)) =
gammasgn(x) * gamma(x)
.
For complex-valued log-gamma, use loggamma
instead of gammaln
.
References
.. [dlmf] NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5
Examples
>>> import numpy as np
>>> import scipy.special as sc
It has two positive zeros.
>>> sc.gammaln([1, 2])
array([0., 0.])
It has poles at nonpositive integers.
>>> sc.gammaln([0, -1, -2, -3, -4])
array([inf, inf, inf, inf, inf])
It asymptotically approaches x * log(x)
(Stirling's formula).
>>> x = np.array([1e10, 1e20, 1e40, 1e80])
>>> sc.gammaln(x)
array([2.20258509e+11, 4.50517019e+21, 9.11034037e+41, 1.83206807e+82])
>>> x * np.log(x)
array([2.30258509e+11, 4.60517019e+21, 9.21034037e+41, 1.84206807e+82])
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
gammainc(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
gammainc(a, x, out=None)
Regularized lower incomplete gamma function.
It is defined as
$$P(a, x) = \frac{1}{\Gamma(a)} \int_0^x t^{a - 1}e^{-t} dt$$
for \( a > 0 \) and \( x \geq 0 \). See [dlmf]_ for details.
Parameters
- a (array_like): Positive parameter
- x (array_like): Nonnegative argument
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Values of the lower incomplete gamma function
See Also
gammaincc
: regularized upper incomplete gamma function
gammaincinv
: inverse of the regularized lower incomplete gamma function
gammainccinv
: inverse of the regularized upper incomplete gamma function
Notes
The function satisfies the relation gammainc(a, x) +
gammaincc(a, x) = 1
where gammaincc
is the regularized upper
incomplete gamma function.
The implementation largely follows that of [boost]_.
References
.. [dlmf] NIST Digital Library of Mathematical functions https://dlmf.nist.gov/8.2#E4 .. [boost] Maddock et. al., "Incomplete Gamma Functions", https://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
Examples
>>> import scipy.special as sc
It is the CDF of the gamma distribution, so it starts at 0 and monotonically increases to 1.
>>> sc.gammainc(0.5, [0, 1, 10, 100])
array([0. , 0.84270079, 0.99999226, 1. ])
It is equal to one minus the upper incomplete gamma function.
>>> a, x = 0.5, 0.4
>>> sc.gammainc(a, x)
0.6289066304773024
>>> 1 - sc.gammaincc(a, x)
0.6289066304773024
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
gammaincc(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
gammaincc(a, x, out=None)
Regularized upper incomplete gamma function.
It is defined as
$$Q(a, x) = \frac{1}{\Gamma(a)} \int_x^\infty t^{a - 1}e^{-t} dt$$
for \( a > 0 \) and \( x \geq 0 \). See [dlmf]_ for details.
Parameters
- a (array_like): Positive parameter
- x (array_like): Nonnegative argument
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Values of the upper incomplete gamma function
See Also
gammainc
: regularized lower incomplete gamma function
gammaincinv
: inverse of the regularized lower incomplete gamma function
gammainccinv
: inverse of the regularized upper incomplete gamma function
Notes
The function satisfies the relation gammainc(a, x) +
gammaincc(a, x) = 1
where gammainc
is the regularized lower
incomplete gamma function.
The implementation largely follows that of [boost]_.
References
.. [dlmf] NIST Digital Library of Mathematical functions https://dlmf.nist.gov/8.2#E4 .. [boost] Maddock et. al., "Incomplete Gamma Functions", https://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
Examples
>>> import scipy.special as sc
It is the survival function of the gamma distribution, so it starts at 1 and monotonically decreases to 0.
>>> sc.gammaincc(0.5, [0, 1, 10, 100, 1000])
array([1.00000000e+00, 1.57299207e-01, 7.74421643e-06, 2.08848758e-45,
0.00000000e+00])
It is equal to one minus the lower incomplete gamma function.
>>> a, x = 0.5, 0.4
>>> sc.gammaincc(a, x)
0.37109336952269756
>>> 1 - sc.gammainc(a, x)
0.37109336952269756
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
gammasgn(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
gammasgn(x, out=None)
Sign of the gamma function.
It is defined as
$$\text{gammasgn}(x) = \begin{cases} +1 & \Gamma(x) > 0 \ -1 & \Gamma(x) < 0 \end{cases}$$
where \( \Gamma \) is the gamma function; see gamma
. This
definition is complete since the gamma function is never zero;
see the discussion after [dlmf]_.
Parameters
- x (array_like): Real argument
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Sign of the gamma function
See Also
gamma
: the gamma function
gammaln
: log of the absolute value of the gamma function
loggamma
: analytic continuation of the log of the gamma function
Notes
The gamma function can be computed as gammasgn(x) *
np.exp(gammaln(x))
.
References
.. [dlmf] NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5.2#E1
Examples
>>> import numpy as np
>>> import scipy.special as sc
It is 1 for x > 0
.
>>> sc.gammasgn([1, 2, 3, 4])
array([1., 1., 1., 1.])
It alternates between -1 and 1 for negative integers.
>>> sc.gammasgn([-0.5, -1.5, -2.5, -3.5])
array([-1., 1., -1., 1.])
It can be used to compute the gamma function.
>>> x = [1.5, 0.5, -0.5, -1.5]
>>> sc.gammasgn(x) * np.exp(sc.gammaln(x))
array([ 0.88622693, 1.77245385, -3.5449077 , 2.3632718 ])
>>> sc.gamma(x)
array([ 0.88622693, 1.77245385, -3.5449077 , 2.3632718 ])
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
rgamma(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
rgamma(z, out=None)
Reciprocal of the gamma function.
Defined as \( 1 / \Gamma(z) \), where \( \Gamma \) is the
gamma function. For more on the gamma function see gamma
.
Parameters
- z (array_like): Real or complex valued input
- out (ndarray, optional): Optional output array for the function results
Returns
- scalar or ndarray: Function results
See Also
gamma,
, gammaln,
, loggamma
Notes
The gamma function has no zeros and has simple poles at
nonpositive integers, so rgamma
is an entire function with zeros
at the nonpositive integers. See the discussion in [dlmf]_ for
more details.
References
.. [dlmf] Nist, Digital Library of Mathematical functions, https://dlmf.nist.gov/5.2#i
Examples
>>> import scipy.special as sc
It is the reciprocal of the gamma function.
>>> sc.rgamma([1, 2, 3, 4])
array([1. , 1. , 0.5 , 0.16666667])
>>> 1 / sc.gamma([1, 2, 3, 4])
array([1. , 1. , 0.5 , 0.16666667])
It is zero at nonpositive integers.
>>> sc.rgamma([0, -1, -2, -3])
array([0., 0., 0., 0.])
It rapidly underflows to zero along the positive real axis.
>>> sc.rgamma([10, 100, 179])
array([2.75573192e-006, 1.07151029e-156, 0.00000000e+000])
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
Returns the log of multivariate gamma, also sometimes called the generalized gamma.
Parameters
- a (ndarray):
The multivariate gamma is computed for each item of
a
. - d (int): The dimension of the space of integration.
Returns
- res (ndarray):
The values of the log multivariate gamma at the given points
a
.
Notes
The formal definition of the multivariate gamma of dimension d for a real
a
is
$$\Gamma_d(a) = \int_{A>0} e^{-tr(A)} |A|^{a - (d+1)/2} dA$$
with the condition \( a > (d-1)/2 \), and \( A > 0 \) being the set of
all the positive definite matrices of dimension d
. Note that a
is a
scalar: the integrand only is multivariate, the argument is not (the
function is defined over a subset of the real set).
This can be proven to be equal to the much friendlier equation
$$\Gamma_d(a) = \pi^{d(d-1)/4} \prod_{i=1}^{d} \Gamma(a - (i-1)/2).$$
References
R. J. Muirhead, Aspects of multivariate statistical theory (Wiley Series in probability and mathematical statistics).
Examples
>>> import numpy as np
>>> from scipy.special import multigammaln, gammaln
>>> a = 23.5
>>> d = 10
>>> multigammaln(a, d)
454.1488605074416
Verify that the result agrees with the logarithm of the equation shown above:
>>> d*(d-1)/4*np.log(np.pi) + gammaln(a - 0.5*np.arange(0, d)).sum()
454.1488605074416
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
Modified Bessel function of the second kind of integer order n
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
j0(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
j0(x, out=None)
Bessel function of the first kind of order 0.
Parameters
- x (array_like): Argument (float).
- out (ndarray, optional): Optional output array for the function values
Returns
- J (scalar or ndarray):
Value of the Bessel function of the first kind of order 0 at
x
.
See Also
jv
: Bessel function of real order and complex argument.
spherical_jn
: spherical Bessel functions.
Notes
The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval the following rational approximation is used:
$$J_0(x) \approx (w - r_1^2)(w - r_2^2) \frac{P_3(w)}{Q_8(w)},$$
where \( w = x^2 \) and \( r_1 \), \( r_2 \) are the zeros of \( J_0 \), and \( P_3 \) and \( Q_8 \) are polynomials of degrees 3 and 8, respectively.
In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7.
This function is a wrapper for the Cephes 1 routine j0
.
It should not be confused with the spherical Bessel functions (see
spherical_jn
).
References
Examples
Calculate the function at one point:
>>> from scipy.special import j0
>>> j0(1.)
0.7651976865579665
Calculate the function at several points:
>>> import numpy as np
>>> j0(np.array([-2., 0., 4.]))
array([ 0.22389078, 1. , -0.39714981])
Plot the function from -20 to 20.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-20., 20., 1000)
>>> y = j0(x)
>>> ax.plot(x, y)
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
y0(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
y0(x, out=None)
Bessel function of the second kind of order 0.
Parameters
- x (array_like): Argument (float).
- out (ndarray, optional): Optional output array for the function results
Returns
- Y (scalar or ndarray):
Value of the Bessel function of the second kind of order 0 at
x
.
See Also
j0
: Bessel function of the first kind of order 0
yv
: Bessel function of the first kind
Notes
The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval a rational approximation \( R(x) \) is employed to compute,
$$Y_0(x) = R(x) + \frac{2 \log(x) J_0(x)}{\pi},$$
where \( J_0 \) is the Bessel function of the first kind of order 0.
In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7.
This function is a wrapper for the Cephes 1 routine y0
.
References
Examples
Calculate the function at one point:
>>> from scipy.special import y0
>>> y0(1.)
0.08825696421567697
Calculate at several points:
>>> import numpy as np
>>> y0(np.array([0.5, 2., 3.]))
array([-0.44451873, 0.51037567, 0.37685001])
Plot the function from 0 to 10.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0., 10., 1000)
>>> y = y0(x)
>>> ax.plot(x, y)
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
j1(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
j1(x, out=None)
Bessel function of the first kind of order 1.
Parameters
- x (array_like): Argument (float).
- out (ndarray, optional): Optional output array for the function values
Returns
- J (scalar or ndarray):
Value of the Bessel function of the first kind of order 1 at
x
.
See Also
jv
: Bessel function of the first kind
spherical_jn
: spherical Bessel functions.
Notes
The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.
This function is a wrapper for the Cephes 1 routine j1
.
It should not be confused with the spherical Bessel functions (see
spherical_jn
).
References
Examples
Calculate the function at one point:
>>> from scipy.special import j1
>>> j1(1.)
0.44005058574493355
Calculate the function at several points:
>>> import numpy as np
>>> j1(np.array([-2., 0., 4.]))
array([-0.57672481, 0. , -0.06604333])
Plot the function from -20 to 20.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-20., 20., 1000)
>>> y = j1(x)
>>> ax.plot(x, y)
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
y1(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
y1(x, out=None)
Bessel function of the second kind of order 1.
Parameters
- x (array_like): Argument (float).
- out (ndarray, optional): Optional output array for the function results
Returns
- Y (scalar or ndarray):
Value of the Bessel function of the second kind of order 1 at
x
.
See Also
j1
: Bessel function of the first kind of order 1
yn
: Bessel function of the second kind
yv
: Bessel function of the second kind
Notes
The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 25 term Chebyshev expansion is used, and computing \( J_1 \) (the Bessel function of the first kind) is required. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.
This function is a wrapper for the Cephes 1 routine y1
.
References
Examples
Calculate the function at one point:
>>> from scipy.special import y1
>>> y1(1.)
-0.7812128213002888
Calculate at several points:
>>> import numpy as np
>>> y1(np.array([0.5, 2., 3.]))
array([-1.47147239, -0.10703243, 0.32467442])
Plot the function from 0 to 10.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0., 10., 1000)
>>> y = y1(x)
>>> ax.plot(x, y)
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
jv(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
jv(v, z, out=None)
Bessel function of the first kind of real order and complex argument.
Parameters
- v (array_like): Order (float).
- z (array_like): Argument (float or complex).
- out (ndarray, optional): Optional output array for the function values
Returns
- J (scalar or ndarray): Value of the Bessel function, \( J_v(z) \).
See Also
jve
: \( J_v \) with leading exponential behavior stripped off.
spherical_jn
: spherical Bessel functions.
j0
: faster version of this function for order 0.
j1
: faster version of this function for order 1.
Notes
For positive v
values, the computation is carried out using the AMOS
1 zbesj
routine, which exploits the connection to the modified
Bessel function \( I_v \),
$$J_v(z) = \exp(v\pi\imath/2) I_v(-\imath z)\qquad (\Im z > 0)
J_v(z) = \exp(-v\pi\imath/2) I_v(\imath z)\qquad (\Im z < 0)$$
For negative v
values the formula,
$$J_{-v}(z) = J_v(z) \cos(\pi v) - Y_v(z) \sin(\pi v)$$
is used, where \( Y_v(z) \) is the Bessel function of the second
kind, computed using the AMOS routine zbesy
. Note that the second
term is exactly zero for integer v
; to improve accuracy the second
term is explicitly omitted for v
values such that v = floor(v)
.
Not to be confused with the spherical Bessel functions (see spherical_jn
).
References
Examples
Evaluate the function of order 0 at one point.
>>> from scipy.special import jv
>>> jv(0, 1.)
0.7651976865579666
Evaluate the function at one point for different orders.
>>> jv(0, 1.), jv(1, 1.), jv(1.5, 1.)
(0.7651976865579666, 0.44005058574493355, 0.24029783912342725)
The evaluation for different orders can be carried out in one call by
providing a list or NumPy array as argument for the v
parameter:
>>> jv([0, 1, 1.5], 1.)
array([0.76519769, 0.44005059, 0.24029784])
Evaluate the function at several points for order 0 by providing an
array for z
.
>>> import numpy as np
>>> points = np.array([-2., 0., 3.])
>>> jv(0, points)
array([ 0.22389078, 1. , -0.26005195])
If z
is an array, the order parameter v
must be broadcastable to
the correct shape if different orders shall be computed in one call.
To calculate the orders 0 and 1 for an 1D array:
>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)
>>> jv(orders, points)
array([[ 0.22389078, 1. , -0.26005195],
[-0.57672481, 0. , 0.33905896]])
Plot the functions of order 0 to 3 from -10 to 10.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> for i in range(4):
... ax.plot(x, jv(i, x), label=f'$J_{i!r}$')
>>> ax.legend()
>>> plt.show()
-
Donald E. Amos, "AMOS, A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order", http://netlib.org/amos/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
yn(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
yn(n, x, out=None)
Bessel function of the second kind of integer order and real argument.
Parameters
- n (array_like): Order (integer).
- x (array_like): Argument (float).
- out (ndarray, optional): Optional output array for the function results
Returns
- Y (scalar or ndarray): Value of the Bessel function, \( Y_n(x) \).
See Also
yv
: For real order and real or complex argument.
y0
: faster implementation of this function for order 0
y1
: faster implementation of this function for order 1
Notes
Wrapper for the Cephes 1 routine yn
.
The function is evaluated by forward recurrence on n
, starting with
values computed by the Cephes routines y0
and y1
. If n = 0
or 1,
the routine for y0
or y1
is called directly.
References
Examples
Evaluate the function of order 0 at one point.
>>> from scipy.special import yn
>>> yn(0, 1.)
0.08825696421567697
Evaluate the function at one point for different orders.
>>> yn(0, 1.), yn(1, 1.), yn(2, 1.)
(0.08825696421567697, -0.7812128213002888, -1.6506826068162546)
The evaluation for different orders can be carried out in one call by
providing a list or NumPy array as argument for the v
parameter:
>>> yn([0, 1, 2], 1.)
array([ 0.08825696, -0.78121282, -1.65068261])
Evaluate the function at several points for order 0 by providing an
array for z
.
>>> import numpy as np
>>> points = np.array([0.5, 3., 8.])
>>> yn(0, points)
array([-0.44451873, 0.37685001, 0.22352149])
If z
is an array, the order parameter v
must be broadcastable to
the correct shape if different orders shall be computed in one call.
To calculate the orders 0 and 1 for an 1D array:
>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)
>>> yn(orders, points)
array([[-0.44451873, 0.37685001, 0.22352149],
[-1.47147239, 0.32467442, -0.15806046]])
Plot the functions of order 0 to 3 from 0 to 10.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0., 10., 1000)
>>> for i in range(4):
... ax.plot(x, yn(i, x), label=f'$Y_{i!r}$')
>>> ax.set_ylim(-3, 1)
>>> ax.legend()
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
i0(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
i0(x, out=None)
Modified Bessel function of order 0.
Defined as,
$$I_0(x) = \sum_{k=0}^\infty \frac{(x^2/4)^k}{(k!)^2} = J_0(\imath x),$$
where \( J_0 \) is the Bessel function of the first kind of order 0.
Parameters
- x (array_like): Argument (float)
- out (ndarray, optional): Optional output array for the function values
Returns
- I (scalar or ndarray):
Value of the modified Bessel function of order 0 at
x
.
See Also
iv
: Modified Bessel function of any order
i0e
: Exponentially scaled modified Bessel function of order 0
Notes
The range is partitioned into the two intervals [0, 8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.
This function is a wrapper for the Cephes 1 routine i0
.
References
Examples
Calculate the function at one point:
>>> from scipy.special import i0
>>> i0(1.)
1.2660658777520082
Calculate at several points:
>>> import numpy as np
>>> i0(np.array([-2., 0., 3.5]))
array([2.2795853 , 1. , 7.37820343])
Plot the function from -10 to 10.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> y = i0(x)
>>> ax.plot(x, y)
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
i1(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
i1(x, out=None)
Modified Bessel function of order 1.
Defined as,
$$I_1(x) = \frac{1}{2}x \sum_{k=0}^\infty \frac{(x^2/4)^k}{k! (k + 1)!} = -\imath J_1(\imath x),$$
where \( J_1 \) is the Bessel function of the first kind of order 1.
Parameters
- x (array_like): Argument (float)
- out (ndarray, optional): Optional output array for the function values
Returns
- I (scalar or ndarray):
Value of the modified Bessel function of order 1 at
x
.
See Also
iv
: Modified Bessel function of the first kind
i1e
: Exponentially scaled modified Bessel function of order 1
Notes
The range is partitioned into the two intervals [0, 8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.
This function is a wrapper for the Cephes 1 routine i1
.
References
Examples
Calculate the function at one point:
>>> from scipy.special import i1
>>> i1(1.)
0.5651591039924851
Calculate the function at several points:
>>> import numpy as np
>>> i1(np.array([-2., 0., 6.]))
array([-1.59063685, 0. , 61.34193678])
Plot the function between -10 and 10.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> y = i1(x)
>>> ax.plot(x, y)
>>> plt.show()
-
Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
iv(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
iv(v, z, out=None)
Modified Bessel function of the first kind of real order.
Parameters
- v (array_like):
Order. If
z
is of real type and negative,v
must be integer valued. - z (array_like of float or complex): Argument.
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Values of the modified Bessel function.
See Also
ive
: This function with leading exponential behavior stripped off.
i0
: Faster version of this function for order 0.
i1
: Faster version of this function for order 1.
Notes
For real z
and \( v \in [-50, 50] \), the evaluation is carried out
using Temme's method 1. For larger orders, uniform asymptotic
expansions are applied.
For complex z
and positive v
, the AMOS 2 zbesi
routine is
called. It uses a power series for small z
, the asymptotic expansion
for large abs(z)
, the Miller algorithm normalized by the Wronskian
and a Neumann series for intermediate magnitudes, and the uniform
asymptotic expansions for \( I_v(z) \) and \( J_v(z) \) for large
orders. Backward recurrence is used to generate sequences or reduce
orders when necessary.
The calculations above are done in the right half plane and continued into the left half plane by the formula,
$$I_v(z \exp(\pm\imath\pi)) = \exp(\pm\pi v) I_v(z)$$
(valid when the real part of z
is positive). For negative v
, the
formula
$$I_{-v}(z) = I_v(z) + \frac{2}{\pi} \sin(\pi v) K_v(z)$$
is used, where \( K_v(z) \) is the modified Bessel function of the
second kind, evaluated using the AMOS routine zbesk
.
References
Examples
Evaluate the function of order 0 at one point.
>>> from scipy.special import iv
>>> iv(0, 1.)
1.2660658777520084
Evaluate the function at one point for different orders.
>>> iv(0, 1.), iv(1, 1.), iv(1.5, 1.)
(1.2660658777520084, 0.565159103992485, 0.2935253263474798)
The evaluation for different orders can be carried out in one call by
providing a list or NumPy array as argument for the v
parameter:
>>> iv([0, 1, 1.5], 1.)
array([1.26606588, 0.5651591 , 0.29352533])
Evaluate the function at several points for order 0 by providing an
array for z
.
>>> import numpy as np
>>> points = np.array([-2., 0., 3.])
>>> iv(0, points)
array([2.2795853 , 1. , 4.88079259])
If z
is an array, the order parameter v
must be broadcastable to
the correct shape if different orders shall be computed in one call.
To calculate the orders 0 and 1 for an 1D array:
>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)
>>> iv(orders, points)
array([[ 2.2795853 , 1. , 4.88079259],
[-1.59063685, 0. , 3.95337022]])
Plot the functions of order 0 to 3 from -5 to 5.
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-5., 5., 1000)
>>> for i in range(4):
... ax.plot(x, iv(i, x), label=f'$I_{i!r}$')
>>> ax.legend()
>>> plt.show()
-
Temme, Journal of Computational Physics, vol 21, 343 (1976) ↩
-
Donald E. Amos, "AMOS, A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order", http://netlib.org/amos/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
ive(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
ive(v, z, out=None)
Exponentially scaled modified Bessel function of the first kind.
Defined as::
ive(v, z) = iv(v, z) * exp(-abs(z.real))
For imaginary numbers without a real part, returns the unscaled
Bessel function of the first kind iv
.
Parameters
- v (array_like of float): Order.
- z (array_like of float or complex): Argument.
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: Values of the exponentially scaled modified Bessel function.
See Also
iv
: Modified Bessel function of the first kind
i0e
: Faster implementation of this function for order 0
i1e
: Faster implementation of this function for order 1
Notes
For positive v
, the AMOS 1 zbesi
routine is called. It uses a
power series for small z
, the asymptotic expansion for large
abs(z)
, the Miller algorithm normalized by the Wronskian and a
Neumann series for intermediate magnitudes, and the uniform asymptotic
expansions for \( I_v(z) \) and \( J_v(z) \) for large orders.
Backward recurrence is used to generate sequences or reduce orders when
necessary.
The calculations above are done in the right half plane and continued into the left half plane by the formula,
$$I_v(z \exp(\pm\imath\pi)) = \exp(\pm\pi v) I_v(z)$$
(valid when the real part of z
is positive). For negative v
, the
formula
$$I_{-v}(z) = I_v(z) + \frac{2}{\pi} \sin(\pi v) K_v(z)$$
is used, where \( K_v(z) \) is the modified Bessel function of the
second kind, evaluated using the AMOS routine zbesk
.
ive
is useful for large arguments z
: for these, iv
easily overflows,
while ive
does not due to the exponential scaling.
References
Examples
In the following example iv
returns infinity whereas ive
still returns
a finite number.
>>> from scipy.special import iv, ive
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> iv(3, 1000.), ive(3, 1000.)
(inf, 0.01256056218254712)
Evaluate the function at one point for different orders by
providing a list or NumPy array as argument for the v
parameter:
>>> ive([0, 1, 1.5], 1.)
array([0.46575961, 0.20791042, 0.10798193])
Evaluate the function at several points for order 0 by providing an
array for z
.
>>> points = np.array([-2., 0., 3.])
>>> ive(0, points)
array([0.30850832, 1. , 0.24300035])
Evaluate the function at several points for different orders by
providing arrays for both v
for z
. Both arrays have to be
broadcastable to the correct shape. To calculate the orders 0, 1
and 2 for a 1D array of points:
>>> ive([[0], [1], [2]], points)
array([[ 0.30850832, 1. , 0.24300035],
[-0.21526929, 0. , 0.19682671],
[ 0.09323903, 0. , 0.11178255]])
Plot the functions of order 0 to 3 from -5 to 5.
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-5., 5., 1000)
>>> for i in range(4):
... ax.plot(x, ive(i, x), label=fr'$I_{i!r}(z)\cdot e^{{-|z|}}$')
>>> ax.legend()
>>> ax.set_xlabel(r"$z$")
>>> plt.show()
-
Donald E. Amos, "AMOS, A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order", http://netlib.org/amos/ ↩
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
erf(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
erf(z, out=None)
Returns the error function of complex argument.
It is defined as 2/sqrt(pi)*integral(exp(-t**2), t=0..z)
.
Parameters
- x (ndarray): Input array.
- out (ndarray, optional): Optional output array for the function values
Returns
- res (scalar or ndarray):
The values of the error function at the given points
x
.
See Also
erfc,
, erfinv,
, erfcinv,
, wofz,
, erfcx,
, erfi
Notes
The cumulative of the unit normal distribution is given by
Phi(z) = 1/2[1 + erf(z/sqrt(2))]
.
References
Examples
>>> import numpy as np
>>> from scipy import special
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-3, 3)
>>> plt.plot(x, special.erf(x))
>>> plt.xlabel('$x$')
>>> plt.ylabel('$erf(x)$')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
erfc(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
erfc(x, out=None)
Complementary error function, 1 - erf(x)
.
Parameters
- x (array_like): Real or complex valued argument
- out (ndarray, optional): Optional output array for the function results
Returns
- scalar or ndarray: Values of the complementary error function
See Also
erf,
, erfi,
, erfcx,
, dawsn,
, wofz
References
Examples
>>> import numpy as np
>>> from scipy import special
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-3, 3)
>>> plt.plot(x, special.erfc(x))
>>> plt.xlabel('$x$')
>>> plt.ylabel('$erfc(x)$')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
erfinv(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
erfinv(y, out=None)
Inverse of the error function.
Computes the inverse of the error function.
In the complex domain, there is no unique complex number w satisfying erf(w)=z. This indicates a true inverse function would be multivalued. When the domain restricts to the real, -1 < x < 1, there is a unique real number satisfying erf(erfinv(x)) = x.
Parameters
- y (ndarray): Argument at which to evaluate. Domain: [-1, 1]
- out (ndarray, optional): Optional output array for the function values
Returns
- erfinv (scalar or ndarray): The inverse of erf of y, element-wise
See Also
erf
: Error function of a complex argument
erfc
: Complementary error function, 1 - erf(x)
erfcinv
: Inverse of the complementary error function
Examples
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.special import erfinv, erf
>>> erfinv(0.5)
0.4769362762044699
>>> y = np.linspace(-1.0, 1.0, num=9)
>>> x = erfinv(y)
>>> x
array([ -inf, -0.81341985, -0.47693628, -0.22531206, 0. ,
0.22531206, 0.47693628, 0.81341985, inf])
Verify that erf(erfinv(y))
is y
.
>>> erf(x)
array([-1. , -0.75, -0.5 , -0.25, 0. , 0.25, 0.5 , 0.75, 1. ])
Plot the function:
>>> y = np.linspace(-1, 1, 200)
>>> fig, ax = plt.subplots()
>>> ax.plot(y, erfinv(y))
>>> ax.grid(True)
>>> ax.set_xlabel('y')
>>> ax.set_title('erfinv(y)')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
erfcinv(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
erfcinv(y, out=None)
Inverse of the complementary error function.
Computes the inverse of the complementary error function.
In the complex domain, there is no unique complex number w satisfying erfc(w)=z. This indicates a true inverse function would be multivalued. When the domain restricts to the real, 0 < x < 2, there is a unique real number satisfying erfc(erfcinv(x)) = erfcinv(erfc(x)).
It is related to inverse of the error function by erfcinv(1-x) = erfinv(x)
Parameters
- y (ndarray): Argument at which to evaluate. Domain: [0, 2]
- out (ndarray, optional): Optional output array for the function values
Returns
- erfcinv (scalar or ndarray): The inverse of erfc of y, element-wise
See Also
erf
: Error function of a complex argument
erfc
: Complementary error function, 1 - erf(x)
erfinv
: Inverse of the error function
Examples
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.special import erfcinv
>>> erfcinv(0.5)
0.4769362762044699
>>> y = np.linspace(0.0, 2.0, num=11)
>>> erfcinv(y)
array([ inf, 0.9061938 , 0.59511608, 0.37080716, 0.17914345,
-0. , -0.17914345, -0.37080716, -0.59511608, -0.9061938 ,
-inf])
Plot the function:
>>> y = np.linspace(0, 2, 200)
>>> fig, ax = plt.subplots()
>>> ax.plot(y, erfcinv(y))
>>> ax.grid(True)
>>> ax.set_xlabel('y')
>>> ax.set_title('erfcinv(y)')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
logit(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
""" logit(x, out=None)
Logit ufunc for ndarrays.
The logit function is defined as logit(p) = log(p/(1-p)). Note that logit(0) = -inf, logit(1) = inf, and logit(p) for p<0 or p>1 yields nan.
Parameters
- x (ndarray): The ndarray to apply logit to element-wise.
- out (ndarray, optional): Optional output array for the function results
Returns
- scalar or ndarray: An ndarray of the same shape as x. Its entries are logit of the corresponding entry of x.
See Also
Notes
As a ufunc logit takes a number of optional keyword arguments. For more information see ufuncs
New in version 0.10.0.
Examples
>>> import numpy as np
>>> from scipy.special import logit, expit
>>> logit([0, 0.25, 0.5, 0.75, 1])
array([ -inf, -1.09861229, 0. , 1.09861229, inf])
expit
is the inverse of logit
:
>>> expit(logit([0.1, 0.75, 0.999]))
array([ 0.1 , 0.75 , 0.999])
Plot logit(x) for x in [0, 1]:
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0, 1, 501)
>>> y = logit(x)
>>> plt.plot(x, y)
>>> plt.grid()
>>> plt.ylim(-6, 6)
>>> plt.xlabel('x')
>>> plt.title('logit(x)')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
expit(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])
expit(x, out=None)
Expit (a.k.a. logistic sigmoid) ufunc for ndarrays.
The expit function, also known as the logistic sigmoid function, is
defined as expit(x) = 1/(1+exp(-x))
. It is the inverse of the
logit function.
Parameters
- x (ndarray): The ndarray to apply expit to element-wise.
- out (ndarray, optional): Optional output array for the function values
Returns
- scalar or ndarray: An ndarray of the same shape as x. Its entries
are
expit
of the corresponding entry of x.
See Also
Notes
As a ufunc expit takes a number of optional keyword arguments. For more information see ufuncs
New in version 0.10.0.
Examples
>>> import numpy as np
>>> from scipy.special import expit, logit
>>> expit([-np.inf, -1.5, 0, 1.5, np.inf])
array([ 0. , 0.18242552, 0.5 , 0.81757448, 1. ])
logit
is the inverse of expit
:
>>> logit(expit([-2.5, 0, 3.1, 5.0]))
array([-2.5, 0. , 3.1, 5. ])
Plot expit(x) for x in [-6, 6]:
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-6, 6, 121)
>>> y = expit(x)
>>> plt.plot(x, y)
>>> plt.grid()
>>> plt.xlim(-6, 6)
>>> plt.xlabel('x')
>>> plt.title('expit(x)')
>>> plt.show()
36 @wraps(f_raw) 37 def f_wrapped(*args, **kwargs): 38 boxed_args, trace, node_constructor = find_top_boxed_args(args) 39 if boxed_args: 40 argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args]) 41 if f_wrapped in notrace_primitives[node_constructor]: 42 return f_wrapped(*argvals, **kwargs) 43 parents = tuple(box._node for _ , box in boxed_args) 44 argnums = tuple(argnum for argnum, _ in boxed_args) 45 ans = f_wrapped(*argvals, **kwargs) 46 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents) 47 return new_box(ans, trace, node) 48 else: 49 return f_raw(*args, **kwargs)
Compute the log of the sum of exponentials of input elements.
Parameters
- a (array_like): Input array.
axis (None or int or tuple of ints, optional): Axis or axes over which the sum is taken. By default
axis
is None, and all elements are summed.New in version 0.11.0.
b (array-like, optional): Scaling factor for exp(
a
) must be of the same shape asa
or broadcastable toa
. These values may be negative in order to implement subtraction.New in version 0.12.0.
keepdims (bool, optional): If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the original array.
New in version 0.15.0.
return_sign (bool, optional): If this is set to True, the result will be a pair containing sign information; if False, results that are negative will be returned as NaN. Default is False (no sign information).
New in version 0.16.0.
Returns
- res (ndarray):
The result,
np.log(np.sum(np.exp(a)))
calculated in a numerically more stable way. Ifb
is given thennp.log(np.sum(b*np.exp(a)))
is returned. Ifreturn_sign
is True,res
contains the log of the absolute value of the argument. - sgn (ndarray):
If
return_sign
is True, this will be an array of floating-point numbers matching res containing +1, 0, -1 (for real-valued inputs) or a complex phase (for complex inputs). This gives the sign of the argument of the logarithm inres
. Ifreturn_sign
is False, only one result is returned.
See Also
numpy.logaddexp,
, numpy.logaddexp2
Notes
NumPy has a logaddexp function which is very similar to logsumexp
, but
only handles two arguments. logaddexp.reduce
is similar to this
function, but may be less stable.
Examples
>>> import numpy as np
>>> from scipy.special import logsumexp
>>> a = np.arange(10)
>>> logsumexp(a)
9.4586297444267107
>>> np.log(np.sum(np.exp(a)))
9.4586297444267107
With weights
>>> a = np.arange(10)
>>> b = np.arange(10, 0, -1)
>>> logsumexp(a, b=b)
9.9170178533034665
>>> np.log(np.sum(b*np.exp(a)))
9.9170178533034647
Returning a sign flag
>>> logsumexp([1,2],b=[1,-1],return_sign=True)
(1.5413248546129181, -1.0)
Notice that logsumexp
does not directly support masked arrays. To use it
on a masked array, convert the mask into zero weights:
>>> a = np.ma.array([np.log(2), 2, np.log(3)],
... mask=[False, True, False])
>>> b = (~a.mask).astype(int)
>>> logsumexp(a.data, b=b), np.log(5)
1.6094379124341005, 1.6094379124341005