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

Code before:

    # all_same_const is True if all the values in the groups along the axis=0
    # slice are the same (e.g. [[3, 3, 3], [3, 3, 3, 3], [3, 3, 3]]).
    all_same_const = xp.all(alldata[..., :1] == alldata, axis=-1)

    if not isinstance(equal_var, bool):
        raise TypeError("Expected a boolean value for 'equal_var'")

    if equal_var:
        # Determine the mean of the data, and subtract that from all inputs to a
        # variance (via sum_of_sq / sq_of_sum) calculation.  Variance is invariant
        # to a shift in location, and centering all data around zero vastly
        # improves numerical stability.
        offset = xp.mean(alldata, axis=-1, keepdims=True)
        alldata = alldata - offset

        normalized_ss = xp.sum(alldata, axis=-1)**2. / bign

        sstot = xp.vecdot(alldata, alldata, axis=-1) - normalized_ss

        ssbn = 0
        for sample in samples:
            smo_ss = xp.sum(sample - offset, axis=-1)**2.
            ssbn = ssbn + smo_ss / sample.shape[-1]

        # Naming: variables ending in bn/b are for "between treatments", wn/w are
        # for "within treatments"
        ssbn = ssbn - normalized_ss
        sswn = sstot - ssbn
        dfbn = num_groups - 1
        dfwn = bign - num_groups
        msb = ssbn / dfbn
        msw = sswn / dfwn
        with np.errstate(divide='ignore', invalid='ignore'):
            f = msb / msw
        dfn, dfd = dfbn, dfwn

    else:
        # calculate basic statistics for each sample
        # Beginning of second paragraph [4] page 1:
        # "As a particular case $y_t$ may be the means ... of samples
        y_t = xp.stack([xp.mean(sample, axis=-1) for sample in samples])
        # "... of $n_t$ observations..."
        n_t = xp.asarray([sample.shape[-1] for sample in samples], dtype=y_t.dtype)
        n_t = xp.reshape(n_t, (-1,) + (1,) * (y_t.ndim - 1))
        # "... from $k$ different normal populations..."
        k = len(samples)
        # "The separate samples provide estimates $s_t^2$ of the $\sigma_t^2$."
        s_t2 = xp.stack([xp.var(sample, axis=-1, correction=True)
                         for sample in samples])


Review comment (suggestion):
Avoiding backend-dependent warnings. Code produces a correct result (`nan` for undefined and `inf` for infinite); no need to add noise.

Code after:
    all_const = xp.all(is_const, axis=-1)

    # all_same_const is True if all the values in the groups along the axis=0
    # slice are the same (e.g. [[3, 3, 3], [3, 3, 3, 3], [3, 3, 3]]).
    all_same_const = xp.all(xp.diff(alldata, axis=-1) == 0, axis=-1)

    if not isinstance(equal_var, bool):
        raise TypeError("Expected a boolean value for 'equal_var'")

    if equal_var:
        # Determine the mean of the data, and subtract that from all inputs to a
        # variance (via sum_of_sq / sq_of_sum) calculation.  Variance is invariant
        # to a shift in location, and centering all data around zero vastly
        # improves numerical stability.
        offset = xp.mean(alldata, axis=-1, keepdims=True)
        alldata = alldata - offset

        normalized_ss = xp.sum(alldata, axis=-1)**2. / bign

        sstot = xp.vecdot(alldata, alldata, axis=-1) - normalized_ss

        ssbn = 0
        for sample in samples:
            smo_ss = xp.sum(sample - offset, axis=-1)**2.
            ssbn = ssbn + smo_ss / _length_nonmasked(sample, axis=-1, xp=xp)

        # Naming: variables ending in bn/b are for "between treatments", wn/w are
        # for "within treatments"
        ssbn = ssbn - normalized_ss
        sswn = sstot - ssbn
        dfbn = num_groups - 1
        dfwn = bign - num_groups
        msb = ssbn / dfbn
        msw = sswn / dfwn
        with np.errstate(divide='ignore', invalid='ignore'):
            f = msb / msw
        dfn, dfd = dfbn, dfwn

    else:
        # calculate basic statistics for each sample
        # Beginning of second paragraph [4] page 1:
        # "As a particular case $y_t$ may be the means ... of samples
        y_t = xp.stack([xp.mean(sample, axis=-1) for sample in samples])
        # "... of $n_t$ observations..."
        if is_marray(xp):
            n_t = xp.stack([_length_nonmasked(sample, axis=-1, xp=xp)
                            for sample in samples])
            n_t = xp.asarray(n_t, dtype=n_t.dtype)
        else:
            n_t = xp.asarray([sample.shape[-1] for sample in samples], dtype=y_t.dtype)
            n_t = xp.reshape(n_t, (-1,) + (1,) * (y_t.ndim - 1))