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

Code before:

    Examples
    --------
    >>> from scipy import stats
    >>> x1 = [1,2,1,1,1,1,1,1,1,1,2,4,1,1]
    >>> x2 = [3,3,4,3,1,2,3,1,1,5,4]
    >>> w, p_value = stats.brunnermunzel(x1, x2)
    >>> w
    3.1374674823029505
    >>> p_value
    0.0057862086661515377

    """
    xp = array_namespace(x, y)
    nx = x.shape[axis]
    ny = y.shape[axis]

    rankc = rankdata(xp.concat((x, y), axis=axis), axis=axis)
    rankcx = rankc[..., 0:nx]
    rankcy = rankc[..., nx:nx+ny]
    rankcx_mean = xp.mean(rankcx, axis=axis, keepdims=True)
    rankcy_mean = xp.mean(rankcy, axis=axis, keepdims=True)
    rankx = rankdata(x, axis=axis)
    ranky = rankdata(y, axis=axis)
    rankx_mean = xp.mean(rankx, axis=axis, keepdims=True)
    ranky_mean = xp.mean(ranky, axis=axis, keepdims=True)

    temp_x = rankcx - rankx - rankcx_mean + rankx_mean
    Sx = xp.vecdot(temp_x, temp_x, axis=-1)
    Sx /= nx - 1
    temp_y = rankcy - ranky - rankcy_mean + ranky_mean
    Sy = xp.vecdot(temp_y, temp_y, axis=-1)
    Sy /= ny - 1

    rankcx_mean = xp.squeeze(rankcx_mean, axis=axis)
    rankcy_mean = xp.squeeze(rankcy_mean, axis=axis)
    wbfn = nx * ny * (rankcy_mean - rankcx_mean)
    wbfn /= (nx + ny) * xp.sqrt(nx * Sx + ny * Sy)

    if distribution == "t":
        df_numer = xp.pow(nx * Sx + ny * Sy, 2.0)
        df_denom = xp.pow(nx * Sx, 2.0) / (nx - 1)
        df_denom += xp.pow(ny * Sy, 2.0) / (ny - 1)
        df = df_numer / df_denom

        if (df_numer == 0) and (df_denom == 0):
            message = ("p-value cannot be estimated with `distribution='t' "
                       "because degrees of freedom parameter is undefined "
                       "(0/0). Try using `distribution='normal'")
            warnings.warn(message, RuntimeWarning, stacklevel=2)


Review comment (suggestion):
```suggestion
    # _axis_nan_policy decorator ensures we can work along the last axis
    rankc = rankdata(xp.concat((x, y), axis=axis), axis=-1)
    rankcx = rankc[..., 0:nx]
    rankcy = rankc[..., nx:nx+ny]
    rankcx_mean = xp.mean(rankcx, axis=-1, keepdims=True)
    rankcy_mean = xp.mean(rankcy, axis=-1, keepdims=True)
    rankx = rankdata(x, axis=-1)
    ranky = rankdata(y, axis=-1)
    rankx_mean = xp.mean(rankx, axis=-1, keepdims=True)
    ranky_mean = xp.mean(ranky, axis=-1, keepdims=True)
```

Code after:
           Biometrical Journal. Vol. 42(2000): 17-25.
    .. [2] Neubert, K. and Brunner, E. "A studentized permutation test for the
           non-parametric Behrens-Fisher problem". Computational Statistics and
           Data Analysis. Vol. 51(2007): 5192-5204.

    Examples
    --------
    >>> from scipy import stats
    >>> x1 = [1,2,1,1,1,1,1,1,1,1,2,4,1,1]
    >>> x2 = [3,3,4,3,1,2,3,1,1,5,4]
    >>> w, p_value = stats.brunnermunzel(x1, x2)
    >>> w
    3.1374674823029505
    >>> p_value
    0.0057862086661515377

    """
    xp = array_namespace(x, y)
    nx = x.shape[-1]
    ny = y.shape[-1]

    # _axis_nan_policy decorator ensures we can work along the last axis
    rankc = rankdata(xp.concat((x, y), axis=axis), axis=-1)
    rankcx = rankc[..., 0:nx]
    rankcy = rankc[..., nx:nx+ny]
    rankcx_mean = xp.mean(rankcx, axis=-1, keepdims=True)
    rankcy_mean = xp.mean(rankcy, axis=-1, keepdims=True)
    rankx = rankdata(x, axis=-1)
    ranky = rankdata(y, axis=-1)
    rankx_mean = xp.mean(rankx, axis=-1, keepdims=True)
    ranky_mean = xp.mean(ranky, axis=-1, keepdims=True)

    temp_x = rankcx - rankx - rankcx_mean + rankx_mean
    Sx = xp.vecdot(temp_x, temp_x, axis=-1)
    Sx /= nx - 1
    temp_y = rankcy - ranky - rankcy_mean + ranky_mean
    Sy = xp.vecdot(temp_y, temp_y, axis=-1)
    Sy /= ny - 1

    rankcx_mean = xp.squeeze(rankcx_mean, axis=-1)
    rankcy_mean = xp.squeeze(rankcy_mean, axis=-1)
    wbfn = nx * ny * (rankcy_mean - rankcx_mean)
    wbfn /= (nx + ny) * xp.sqrt(nx * Sx + ny * Sy)

    if distribution == "t":
        df_numer = xp.pow(nx * Sx + ny * Sy, 2.0)
        df_denom = xp.pow(nx * Sx, 2.0) / (nx - 1)
        df_denom += xp.pow(ny * Sy, 2.0) / (ny - 1)
        df = df_numer / df_denom

        if xp.any(df_numer == 0) and xp.any(df_denom == 0):