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

Code before:
    tn = 1 / 6 + (t - et) / math.sqrt(45 * vt)

    # approximate distribution of tn with limiting distribution
    # of the one-sample test statistic
    # if tn < 0.003, the _cdf_cvm_inf(tn) < 1.28*1e-18, return 1.0 directly
    p = xpx.apply_where(tn >= 0.003,
                        (tn,),
                        lambda tn: xp.clip(1. - _cdf_cvm_inf(tn, xp=xp), 0.),
                        fill_value = 1.)
    return p


@xp_capabilities(skip_backends=[('cupy', 'needs rankdata'),
                                ('dask.array', 'needs rankdata')],
                 cpu_only=True, jax_jit=False)
@_axis_nan_policy_factory(CramerVonMisesResult, n_samples=2, too_small=1,
                          result_to_tuple=_cvm_result_to_tuple)
def cramervonmises_2samp(x, y, method='auto', *, axis=0):
    r"""Perform the two-sample Cramér-von Mises test for goodness of fit.

    This is the two-sample version of the Cramér-von Mises test ([1]_):
    for two independent samples :math:`X_1, ..., X_n` and
    :math:`Y_1, ..., Y_m`, the null hypothesis is that the samples
    come from the same (unspecified) continuous distribution.

    The statistic :math:`T` is defined as in [2]_:

    .. math::
        T = \frac{nm}{n+m}\omega^2 =
        \frac{U}{n m (n+m)} - \frac{4 m n - 1}{6(m+n)}

    where :math:`U` is defined as below. The function :math:`r()` here denotes
    the rank of the element within the pooled sample of size :math:`n + m`, with
    ties assigned mid-rank values:

    .. math::
        U = n \sum_{i=1}^n (r(X_i)-i)^2 + m \sum_{j=1}^m (r(Y_j)-j)^2

    Parameters
    ----------
    x : array_like
        A 1-D array of observed values of the random variables :math:`X_i`.
        Must contain at least two observations.
    y : array_like
        A 1-D array of observed values of the random variables :math:`Y_i`.
        Must contain at least two observations.
    method : {'auto', 'asymptotic', 'exact'}, optional
        The method used to compute the p-value, see Notes for details.
        The default is 'auto'.
    axis : int or tuple of ints, default: 0
        If an int or tuple of ints, the axis or axes of the input along which

Review comment (suggestion):
I don't see the statistic labeled as $T$ in [2]. [2] does define a $U$, but it is quite different from what is here, which would be confusing. Since this is actually what is on Wikipedia, which is more accessible than [2], how about:
```suggestion
    The statistic :math:`T` is defined as in [1]_:
```

Code after:
        for v, g in enumerate(gs):
            # Calculate g recursively with eq. 11 in [1]. Even though it
            # doesn't look like it, this also does 12/13 (all of Algorithm 1).
            vi, i0, i1 = np.intersect1d(tmp[0], g[0], return_indices=True)
            tmp = np.concatenate([
                np.stack([vi, tmp[1, i0] + g[1, i1]]),
                np.delete(tmp, i0, 1),
                np.delete(g, i1, 1)
            ], 1)
            res = (a * v - b * u) ** 2
            tmp[0] += res.astype(dtype)
            next_gs.append(tmp)
        gs = next_gs
    value, freq = gs[m]
    return np.float64(np.sum(freq[value >= zeta]) / combinations)


def _pval_cvm_2samp_asymptotic(t, N, nx, ny, k, *, xp):
    # compute expected value and variance of T (eq. 11 and 14 in [2])
    et = (1 + 1 / N) / 6
    vt = (N + 1) * (4 * k * N - 3 * (nx ** 2 + ny ** 2) - 2 * k)
    vt = vt / (45 * N ** 2 * 4 * k)

    # computed the normalized statistic (eq. 15 in [2])
    tn = 1 / 6 + (t - et) / math.sqrt(45 * vt)

    # approximate distribution of tn with limiting distribution
    # of the one-sample test statistic
    # if tn < 0.003, the _cdf_cvm_inf(tn) < 1.28*1e-18, return 1.0 directly
    p = xpx.apply_where(tn >= 0.003,
                        (tn,),
                        lambda tn: xp.clip(1. - _cdf_cvm_inf(tn, xp=xp), 0.),
                        fill_value = 1.)
    return p


@xp_capabilities(skip_backends=[('cupy', 'needs rankdata'),
                                ('dask.array', 'needs rankdata')],
                 cpu_only=True, jax_jit=False)
@_axis_nan_policy_factory(CramerVonMisesResult, n_samples=2, too_small=1,
                          result_to_tuple=_cvm_result_to_tuple)
def cramervonmises_2samp(x, y, method='auto', *, axis=0):
    r"""Perform the two-sample Cramér-von Mises test for goodness of fit.

    This is the two-sample version of the Cramér-von Mises test ([1]_):
    for two independent samples :math:`X_1, ..., X_n` and
    :math:`Y_1, ..., Y_m`, the null hypothesis is that the samples
    come from the same (unspecified) continuous distribution.

    The test statistic :math:`T` is defined as in [1]_: