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

Code before:
    be better to use permuted Brunner Munzel test (see [2]_).

    References
    ----------
    .. [1] Brunner, E. and Munzel, U. "The nonparametric Benhrens-Fisher
           problem: Asymptotic theory and a small-sample approximation".
           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[axis]
    ny = y.shape[axis]

    # _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":

Review comment (suggestion):
```suggestion
    nx = x.shape[-1]
    ny = y.shape[-1]
```

Code after:
    Notes
    -----
    Brunner and Munzel recommended to estimate the p-value by t-distribution
    when the size of data is 50 or less. If the size is lower than 10, it would
    be better to use permuted Brunner Munzel test (see [2]_).

    References
    ----------
    .. [1] Brunner, E. and Munzel, U. "The nonparametric Benhrens-Fisher
           problem: Asymptotic theory and a small-sample approximation".
           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)