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

Code before:
    allZij = []
    g = [0]
    for i in range(k):
        allZij.extend(list(Zij[i]))
        g.append(len(allZij))

    ranks = _stats_py.rankdata(allZij)
    sample = distributions.norm.ppf(ranks / (2*(Ntot + 1.0)) + 0.5)

    # compute Aibar
    Aibar = _apply_func(sample, g, np.sum) / Ni
    anbar = np.mean(sample, axis=0)
    varsq = np.var(sample, axis=0, ddof=1)
    statistic = np.sum(Ni * (asarray(Aibar) - anbar)**2.0, axis=0) / varsq
    chi2 = _SimpleChi2(k-1)
    pval = _get_pvalue(statistic, chi2, alternative='greater', symmetric=False, xp=np)
    return FlignerResult(statistic, pval)


def _mood_statistic_with_ties(x, y, t, n, m, N):
    # First equation of "Mood's Squared Rank Test", Mielke pg 313
    E_0_T = n * (N * N - 1) / 12

    # m, n, N, t, and S are defined in the second paragraph of Mielke pg 312
    # The only difference is that our `t` has zeros interspersed with the relevant
    # numbers to keep the array rectangular, but these terms nothing to the sum.
    np = array_namespace(x, y, t)  # needed to use `cumulative_sum`
    S = np.cumulative_sum(t, include_initial=True, axis=-1)
    S_i, S_i_m1 = S[..., 1:], S[..., :-1]
    # Second equation of "Mood's Squared Rank Test", Mielke pg 313
    varM = (m * n * (N + 1.0) * (N**2 - 4) / 180
            - m * n / (180 * N * (N - 1))
            * np.sum(t * (t ** 2 - 1) * (t ** 2 - 4 + (15 * (N - S_i - S_i_m1) ** 2)),
                     axis=-1))

    # There is a formula for phi in terms of t, S, and Psi(I) = [I - (N+1)/2]^2 is
    # defined (with a mistake in the location of the ^2) at the beginning of
    # "Mood's Squared Rank Test", Mielke pg 313. To vectorize this calculation,
    # let c = (N + 1) / 2, so Psi(I) = I^2 - 2*c*I + c^2. We sum each of these three
    # parts of Psi separately using formulas for sums from a to b (inclusive) of
    # terms I^2, I, and 1 where I takes on successive integers.
    def sum_I2(a, b):
        return (b * (b + 1) * (2*b + 1) - a * (a - 1) * (2*a - 1)) / 6

    def sum_I(a, b):
        return (b * (b + 1) - a * (a - 1)) / 2

    def sum_1(a, b):
        return (b - a) + 1

    with np.errstate(invalid='ignore', divide='ignore'):

Review comment (suggestion):
```suggestion
    # numbers to keep the array rectangular, but these terms do not change the sum.
```

Code after:
    >>> c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
    >>> stat, p = stats.fligner(a, b, c)
    >>> p
    0.00450826080004775

    The small p-value suggests that the populations do not have equal
    variances.

    This is not surprising, given that the sample variance of `b` is much
    larger than that of `a` and `c`:

    >>> [np.var(x, ddof=1) for x in [a, b, c]]
    [0.007054444444444413, 0.13073888888888888, 0.008890000000000002]

    For a more detailed example, see :ref:`hypothesis_fligner`.
    """
    xp = array_namespace(*samples)

    if center not in ['mean', 'median', 'trimmed']:
        raise ValueError("center must be 'mean', 'median' or 'trimmed'.")

    k = len(samples)
    if k < 2:
        raise ValueError("Must provide at least two samples.")

    samples = xp_promote(*samples, force_floating=True, xp=xp)
    dtype = samples[0].dtype

    # Handle empty input
    for sample in samples:
        if sample.size == 0:
            NaN = _get_nan(*samples, xp=xp)
            return FlignerResult(NaN, NaN)

    if center == 'median':

        def func(x):
            return (xp.median(x, axis=-1, keepdims=True)
                    if (is_numpy(xp) or is_dask(xp))
                    else stats.quantile(x, 0.5, axis=-1, keepdims=True))

    elif center == 'mean':

        def func(x):
            return xp.mean(x, axis=-1, keepdims=True)

    else:  # center == 'trimmed'

        def func(x):
            # keepdims=True doesn't currently work for lazy arrays
            return _stats_py.trim_mean(x, proportiontocut, axis=-1)[..., xp.newaxis]