Repository: scipy/scipy
File: scipy/stats/tests/test_stats.py (Python)

Code before:
        res = stats.quantile_test(x, q=q, p=p, alternative=alternative, axis=axis)
        res_statistic, res_pvalue = res.statistic, res.pvalue
        res_low, res_high = res.confidence_interval()
        if axis == 0:
            res_statistic, res_pvalue = res_statistic.T, res_pvalue.T
            res_low, res_high = res_low.T, res_high.T

        assert_allclose(res_statistic, ref_statistic)
        assert_allclose(res_pvalue, ref_pvalue)
        assert_allclose(res_low, res_low)
        assert_allclose(res_low, ref_low)

    def test_zero_size(self):
        rng = np.random.default_rng(883771738488451943)
        x_shape = (2, 50)
        qp_shape = (2, 10)
        x = rng.random(x_shape)
        q = rng.random(qp_shape)
        p = q + rng.random(qp_shape) * 1e-2

        # case 1: p/q is size zero.
        qp_zero = np.empty((0, *qp_shape))
        out = qp_zero
        res = stats.quantile_test(x, q=qp_zero, p=qp_zero, axis=-1)
        ci = res.confidence_interval()
        xp_assert_equal(res.statistic, np.astype(out, np.int64))
        xp_assert_equal(res.statistic_type, np.astype(out, np.int64))
        xp_assert_equal(res.pvalue, out)
        xp_assert_equal(ci.low, out)
        xp_assert_equal(ci.high, out)

        # case 2: x is size zero with nonzero length along axis.
        x_zero = np.empty((0, *x_shape))
        out = np.empty((0, *qp_shape))
        res = stats.quantile_test(x_zero, q=q, p=p, axis=-1)
        ci = res.confidence_interval()
        xp_assert_equal(res.statistic, np.astype(out, np.int64))
        xp_assert_equal(res.statistic_type, np.astype(out, np.int64))
        xp_assert_equal(res.pvalue, out)
        xp_assert_equal(ci.low, out)
        xp_assert_equal(ci.high, out)

        # case 3: x has zero length along axis.
        x_zero = np.empty((x.shape[0], 0))
        out = -np.ones(qp_shape, dtype=np.int64)
        res = stats.quantile_test(x_zero, q=q, p=p, axis=-1)
        ci = res.confidence_interval()
        xp_assert_equal(res.statistic, out)
        xp_assert_equal(res.statistic_type, out)
        xp_assert_equal(res.pvalue, out*np.nan)
        xp_assert_equal(ci.low, out*np.nan)

Review comment (suggestion):
I used `assert_close`/`assert_equal` elsewhere... no harm done in accidentally using `xp_assert_`s here because that will be the next step anyway.

Code after:
        X = xp.asarray([1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1, xp.nan])
        Y = xp.asarray([3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4.])
        alternative = "two-sided"
        distribution = "t"
        nan_policy = "raise"

        message = "The input contains nan values"
        with pytest.raises(ValueError, match=message):
            stats.brunnermunzel(X, Y, alternative, distribution, nan_policy)
        with pytest.raises(ValueError, match=message):
            stats.brunnermunzel(Y, X, alternative, distribution, nan_policy)

    @pytest.mark.skip_xp_backends('jax.numpy', reason='lazy -> no nan_policy')
    def test_brunnermunzel_nan_input_omit(self, xp):
        X = xp.asarray([1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1, np.nan])
        Y = xp.asarray([3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4.])
        u1, p1 = stats.brunnermunzel(X, Y, nan_policy="omit")
        u2, p2 = stats.brunnermunzel(Y, X, nan_policy="omit")

        xp_assert_close(p1, p2)
        xp_assert_close(u1, xp.asarray(3.1374674823029505))
        xp_assert_close(u2, xp.asarray(-3.1374674823029505))
        xp_assert_close(p1, xp.asarray(0.0057862086661515377))

    @pytest.mark.parametrize("vectorized_call", [False, True])
    def test_brunnermunzel_return_nan(self, vectorized_call, xp):
        """ tests that a warning is emitted when p is nan
        p-value with t-distributions can be nan (0/0) (see gh-15843)
        """
        x = xp.asarray([1., 2., 3.])
        y = xp.asarray([5., 6., 7., 8., 9.])

        if vectorized_call:
            x = xp.stack((x, x)).T
            y = xp.stack((y, y)).T

        msg = "p-value cannot be estimated|divide by zero|invalid value encountered"
        with eager_warns(RuntimeWarning, match=msg, xp=xp):
            stats.brunnermunzel(x, y, distribution="t")

    def test_brunnermunzel_normal_dist(self, xp):
        """ tests that a p is 0 for datasets that cause p->nan
        when t-distribution is used (see gh-15843)
        """
        x = xp.asarray([1., 2., 3.])
        y = xp.asarray([5., 6., 7., 8., 9.])

        with eager_warns(RuntimeWarning, match='divide by zero', xp=xp):
            _, p = stats.brunnermunzel(x, y, distribution="normal")
        xp_assert_equal(p, xp.asarray(0.))