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

Code before:
    data = xp.asarray(data, dtype=dtype)

    n = data.shape[axis]
    if n == 0:
        return _get_nan(data, xp=xp)
    eps = xp.finfo(dtype).eps
    pos = data >= 0  # binary mask

    # There exists numerical instability when abs(lmb) or abs(lmb - 2) is very small
    if xp.all(pos):
        if abs(lmb) < eps:
            logvar = xp.log(xp.var(xp.log1p(data), axis=axis))
        else:
            logvar = _log_var(lmb * xp.log1p(data), xp, axis) - 2 * math.log(abs(lmb))

    elif xp.all(~pos):
        if abs(lmb - 2) < eps:
            logvar = xp.log(xp.var(xp.log1p(-data), axis=axis))
        else:
            logvar = _log_var((2 - lmb) * xp.log1p(-data), xp, axis) - 2 * math.log(
                abs(2 - lmb)
            )

    else:  # mixed positive and negative data
        y = _yeojohnson_transform(data, lmb, xp=xp)
        sigma = xp.var(y, axis=axis)

        # Suppress RuntimeWarning raised by np.log when the variance is too low
        finite_variance = sigma >= xp.finfo(sigma.dtype).smallest_normal
        logvar = xpx.apply_where(finite_variance, (sigma,), xp.log, fill_value=-xp.inf)

    loglike = (-n / 2 * logvar
               + (lmb - 1) * xp.sum(xp.sign(data) * xp.log1p(xp.abs(data)), axis=axis))

    return loglike


@xp_capabilities(np_only=True)
def yeojohnson_normmax(x, brack=None):
    """Compute optimal Yeo-Johnson transform parameter.

    Compute optimal Yeo-Johnson transform parameter for input data, using
    maximum likelihood estimation.

    Parameters
    ----------
    x : array_like
        Input array.
    brack : 2-tuple, optional
        The starting interval for a downhill bracket search with
        `optimize.brent`. Note that this is in most cases not critical; the

Review comment (bug):
```suggestion
    # overflow/precision loss not reported for mixed data; calculate `logvar` directly
    else:  # mixed positive and negative data
```

Code after:
    data = xp.asarray(data, dtype=dtype)

    n = data.shape[axis]
    if n == 0:
        return _get_nan(data, xp=xp)
    eps = xp.finfo(dtype).eps
    # Special case all-positive/negative data to avoid overflow and precision loss
    pos = data >= 0  # binary mask

    # There exists numerical instability when abs(lmb) or abs(lmb - 2) is very small
    if xp.all(pos):
        if abs(lmb) < eps:
            logvar = xp.log(xp.var(xp.log1p(data), axis=axis))
        else:
            logvar = _log_var(lmb * xp.log1p(data), xp, axis) - 2 * math.log(abs(lmb))

    elif xp.all(~pos):
        if abs(lmb - 2) < eps:
            logvar = xp.log(xp.var(xp.log1p(-data), axis=axis))
        else:
            logvar = _log_var((2 - lmb) * xp.log1p(-data), xp, axis) - 2 * math.log(
                abs(2 - lmb)
            )

    # overflow/precision loss not reported for mixed data; calculate `logvar` directly
    else:  # mixed positive and negative data
        y = _yeojohnson_transform(data, lmb, xp=xp)
        sigma = xp.var(y, axis=axis)

        # Suppress RuntimeWarning raised by np.log when the variance is too low
        finite_variance = sigma >= xp.finfo(sigma.dtype).smallest_normal
        logvar = xpx.apply_where(finite_variance, (sigma,), xp.log, fill_value=-xp.inf)

    loglike = (-n / 2 * logvar
               + (lmb - 1) * xp.sum(xp.sign(data) * xp.log1p(xp.abs(data)), axis=axis))

    return loglike


@xp_capabilities(np_only=True)
def yeojohnson_normmax(x, brack=None):
    """Compute optimal Yeo-Johnson transform parameter.

    Compute optimal Yeo-Johnson transform parameter for input data, using
    maximum likelihood estimation.

    Parameters
    ----------
    x : array_like
        Input array.
    brack : 2-tuple, optional