# pyerrors.special

View Source

     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)))

@wraps(f_raw)
def beta(*args, **kwargs): View Source

    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

1. NIST Digital Library of Mathematical Functions, Eq. 5.12.1. https://dlmf.nist.gov/5.12 ↩

@wraps(f_raw)
def betainc(*args, **kwargs): View Source

    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

1. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/8.17 ↩

@wraps(f_raw)
def betaln(*args, **kwargs): View Source

    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

@wraps(f_raw)
def polygamma(*args, **kwargs): View Source

    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

`digamma`

###### 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)

@wraps(f_raw)
def psi(*args, **kwargs): View Source

    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)

1. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩

2. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩

3. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩

4. Fredrik Johansson and others. "mpmath: a Python library for arbitrary-precision floating-point arithmetic" (Version 0.19) http://mpmath.org/ ↩

@wraps(f_raw)
def digamma(*args, **kwargs): View Source

    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)

1. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩

2. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩

3. NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5 ↩

4. Fredrik Johansson and others. "mpmath: a Python library for arbitrary-precision floating-point arithmetic" (Version 0.19) http://mpmath.org/ ↩

@wraps(f_raw)
def gamma(*args, **kwargs): View Source

    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()

@wraps(f_raw)
def gammaln(*args, **kwargs): View Source

    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])

@wraps(f_raw)
def gammainc(*args, **kwargs): View Source

    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

@wraps(f_raw)
def gammaincc(*args, **kwargs): View Source

    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

@wraps(f_raw)
def gammasgn(*args, **kwargs): View Source

    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 ])

@wraps(f_raw)
def rgamma(*args, **kwargs): View Source

    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])

@wraps(f_raw)
def multigammaln(*args, **kwargs): View Source

    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

@wraps(f_raw)
def kn(*args, **kwargs): View Source

    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

@wraps(f_raw)
def j0(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def y0(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def j1(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def y1(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def jn(*args, **kwargs): View Source

    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()

1. Donald E. Amos, "AMOS, A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order", http://netlib.org/amos/ ↩

@wraps(f_raw)
def yn(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def i0(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def i1(*args, **kwargs): View Source

    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()

1. Cephes Mathematical Functions Library, http://www.netlib.org/cephes/ ↩

@wraps(f_raw)
def iv(*args, **kwargs): View Source

    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()

1. Temme, Journal of Computational Physics, vol 21, 343 (1976) ↩

2. Donald E. Amos, "AMOS, A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order", http://netlib.org/amos/ ↩

@wraps(f_raw)
def ive(*args, **kwargs): View Source

    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()

1. Donald E. Amos, "AMOS, A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order", http://netlib.org/amos/ ↩

@wraps(f_raw)
def erf(*args, **kwargs): View Source

    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()

@wraps(f_raw)
def erfc(*args, **kwargs): View Source

    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()

@wraps(f_raw)
def erfinv(*args, **kwargs): View Source

    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()

@wraps(f_raw)
def erfcinv(*args, **kwargs): View Source

    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()

@wraps(f_raw)
def logit(*args, **kwargs): View Source

    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

`expit`

###### 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()

@wraps(f_raw)
def expit(*args, **kwargs): View Source

    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

`logit`

###### 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()

@wraps(f_raw)
def logsumexp(*args, **kwargs): View Source

    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 as `a` or broadcastable to `a`. 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. If `b` is given then `np.log(np.sum(b*np.exp(a)))` is returned. If `return_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 in `res`. If `return_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