Repository: scipy/scipy
File: scipy/stats/_continuous_distns.py (Python)

Code before:
        if p < 1:
            return b / (np.sqrt((p - 1)**2 + b**2) + 1 - p)
        else:
            return (np.sqrt((1 - p)**2 + b**2) - (1 - p)) / b

    def _munp(self, n, p, b):
        num = sc.kve(p + n, b)
        denom = sc.kve(p, b)
        inf_vals = np.isinf(num) | np.isinf(denom)
        if inf_vals.any():
            msg = ("Infinite values encountered in the moment calculation "
                   "involving scipy.special.kve. Values replaced by NaN to "
                   "avoid incorrect results.")
            warnings.warn(msg, RuntimeWarning, stacklevel=3)
            m = np.full_like(num, np.nan, dtype=np.float64)
            m[~inf_vals] = num[~inf_vals] / denom[~inf_vals]
        else:
            m = num / denom
        return m


geninvgauss = geninvgauss_gen(a=0.0, name="geninvgauss")


@np.vectorize(otypes=[np.float64])
def _norminvgauss_quadrature(x, y, a, b):
    mean = b / np.sqrt(a**2 - b**2)
    if np.isneginf(x) and y > abs(mean):
        return 1 - integrate.quad(norminvgauss._pdf, y, np.inf, args=(a, b))[0]
    if np.isposinf(y) and x < -abs(mean):
        return 1 - integrate.quad(norminvgauss._pdf, -np.inf, x, args=(a, b))[0]
    else:
        res = integrate.quad(norminvgauss._pdf, x, y, args=(a, b))[0]
    return np.clip(res, 0, 1)


class norminvgauss_gen(rv_continuous):
    r"""A Normal Inverse Gaussian continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `norminvgauss` is:

    .. math::

        f(x, a, b) = \frac{a \, K_1(a \sqrt{1 + x^2})}{\pi \sqrt{1 + x^2}} \,
                     \exp(\sqrt{a^2 - b^2} + b x)

    where :math:`x` is a real number, the parameter :math:`a` is the tail

Review comment (suggestion):
Maybe add the below to comment on the function

```python
"""
norminvgauss CDF difference over [x, y] via quadrature.
Uses the mean as a tail-switching threshold to preserve precision.
"""
```

Code after:
        if p < 1:
            return b / (np.sqrt((p - 1)**2 + b**2) + 1 - p)
        else:
            return (np.sqrt((1 - p)**2 + b**2) - (1 - p)) / b

    def _munp(self, n, p, b):
        num = sc.kve(p + n, b)
        denom = sc.kve(p, b)
        inf_vals = np.isinf(num) | np.isinf(denom)
        if inf_vals.any():
            msg = ("Infinite values encountered in the moment calculation "
                   "involving scipy.special.kve. Values replaced by NaN to "
                   "avoid incorrect results.")
            warnings.warn(msg, RuntimeWarning, stacklevel=3)
            m = np.full_like(num, np.nan, dtype=np.float64)
            m[~inf_vals] = num[~inf_vals] / denom[~inf_vals]
        else:
            m = num / denom
        return m


geninvgauss = geninvgauss_gen(a=0.0, name="geninvgauss")


@np.vectorize(otypes=[np.float64])
def _norminvgauss_quadrature(x, y, a, b):
    # gh-23196 reported that the norminvgauss CDF would drop to zero in the far right
    # tail. The SF had a similar problem, dropping to zero at the far left.
    # This fixes the bug by using 1 - SF to compute the CDF in the right tail and
    # 1 - CDF to compute the SF in the left tail. The mean is guaranteed to be beyond
    # the median, so there is no loss of precision due to subtractive cancellation.
    mean = b / np.sqrt((a + b) * (a - b))
    if np.isneginf(x) and y > abs(mean):
        return 1 - integrate.quad(norminvgauss._pdf, y, np.inf, args=(a, b))[0]
    if np.isposinf(y) and x < -abs(mean):
        return 1 - integrate.quad(norminvgauss._pdf, -np.inf, x, args=(a, b))[0]
    else:
        res = integrate.quad(norminvgauss._pdf, x, y, args=(a, b))[0]
    return np.clip(res, 0, 1)


class norminvgauss_gen(rv_continuous):
    r"""A Normal Inverse Gaussian continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `norminvgauss` is:

    .. math::