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

Code before:
        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
        dtypemix = xp.result_type(dtype, 1j)
        logyj = xp.zeros_like(data, dtype=dtypemix)

        # x >= 0
        if abs(lmb) < eps:
            logyj = xpx.at(logyj)[pos].set(xp.log(xp.log1p(data[pos]) + 0j))
        else:  # lmbda != 0
            logm1_pos = xp.full_like(data[pos], pi * 1j, dtype=dtypemix)
            logyj = xpx.at(logyj)[pos].set(
                special.logsumexp([lmb * xp.log1p(data[pos]), logm1_pos], axis=0)
                - xp.log(lmb + 0j)
            )

        # x < 0
        if abs(lmb - 2) < eps:
            logyj = xpx.at(logyj)[~pos].set(xp.log(-xp.log1p(-data[~pos]) + 0j))
        else:  # lmbda != 2
            logm1_neg = xp.full_like(data[~pos], pi * 1j, dtype=dtypemix)
            logyj = xpx.at(logyj)[~pos].set(
                special.logsumexp(
                    [(2 - lmb) * xp.log1p(-data[~pos]), logm1_neg], axis=0
                )
                - xp.log(lmb - 2 + 0j)
            )

        logvar = _log_var(logyj, xp, axis)

    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

Review comment (suggestion):
I used to do [`logsumexp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.logsumexp.html) subtraction using $\pi j$, but that was because the `b` parameter (weights) wasn't behaving nicely for vectorized code. I think that's fixed now. If you are going to have to work with complex numbers anyway here, then this is fine but if this is the only reason you need to go complex, consider staying real with negative `b`. If that gets much more complicated, never mind.

Code after:

    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
        The starting interval for a downhill bracket search with
        `optimize.brent`. Note that this is in most cases not critical; the
        final result is allowed to be outside this bracket. If None,
        `optimize.fminbound` is used with bounds that avoid overflow.

    Returns
    -------
    maxlog : float
        The optimal transform parameter found.

    See Also
    --------
    yeojohnson, yeojohnson_llf, yeojohnson_normplot

    Notes