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

Code before:
         (0, "gmean", 0),
         (1, "hmean", -1)]
    )
    def test_weights_2d(self, axis, fun_name, p, xp):
        if fun_name == 'wpmean_reference':
            def fun(a, axis, weights):
                return TestPMean.wpmean_reference(a, p, weights)
        else:
            fun = getattr(stats, fun_name)
        a = np.array([[2, 5], [10, 5], [6, 5]])
        weights = np.array([[10, 1], [5, 1], [3, 1]])
        desired = fun(a, axis=axis, weights=weights)
        check_equal_pmean(a, p, desired, axis=axis, weights=weights, rtol=1e-5, xp=xp)

    def test_infinite_p_gh23111(self):
        # gh-23111 reported that `pmean` didn't work properly with infinite `p`;
        # check that this raises an appropriate error message
        message = "Power mean only implemented for finite `p`"
        with pytest.raises(NotImplementedError, match=message):
            stats.pmean([2], np.inf)


    @pytest.mark.parametrize(
        "p, weights, ref",
        [
            (2e-16, None, 47.23984586445855),
            (2e-16, [1., 2., 1., 1., 1., 2., 1.], 62.11911811020747),
        ],
    )

    def test_small_p_gh23407(self, xp, p, weights, ref):
        # gh-23407 reported that pmean had poor numerical behavior for very small nonzero p
        x = xp.asarray([50., 100., 300., 70., 25., 100., 2.])
        w = None if weights is None else xp.asarray(weights)

        res = stats.pmean(x, p, weights=w)
        xp_assert_close(res, xp.asarray(ref), rtol=1e-14)


@make_xp_test_case(stats.gstd)
class TestGSTD:
    # must add 1 as `gstd` is only defined for positive values
    array_1d = (np.arange(2 * 3 * 4) + 1).tolist()
    gstd_array_1d = 2.294407613602
    array_3d = np.reshape(array_1d, (2, 3, 4)).tolist()

    def test_1d_array(self, xp):
        gstd_actual = stats.gstd(xp.asarray(self.array_1d))
        xp_assert_close(gstd_actual, xp.asarray(self.gstd_array_1d))

    @skip_xp_backends(np_only=True, reason="Only NumPy supports array-like input")

Review comment (suggestion):
```suggestion
            (2e-16, None, 47.23984586445856),
```

I'd prefer to use a reference value from code that is very easy to verify.
```python3
import numpy as np
from mpmath import mp
mp.dps = 1000
x = [50., 100., 300., 70., 25., 100., 2.]
x = np.asarray([mp.mpf(x_) for x_ in x])
p = mp.mpf(2e-16)
res = np.average(x**p)**(1/p)
float(res)  # 47.23984586445856
```
Please add this code as a comment in the test so there is a record of where the reference values came from.

Code after:
         (0, "gmean", 0),
         (1, "hmean", -1)]
    )
    def test_weights_2d(self, axis, fun_name, p, xp):
        if fun_name == 'wpmean_reference':
            def fun(a, axis, weights):
                return TestPMean.wpmean_reference(a, p, weights)
        else:
            fun = getattr(stats, fun_name)
        a = np.array([[2, 5], [10, 5], [6, 5]])
        weights = np.array([[10, 1], [5, 1], [3, 1]])
        desired = fun(a, axis=axis, weights=weights)
        check_equal_pmean(a, p, desired, axis=axis, weights=weights, rtol=1e-5, xp=xp)

    def test_infinite_p_gh23111(self):
        # gh-23111 reported that `pmean` didn't work properly with infinite `p`;
        # check that this raises an appropriate error message
        message = "Power mean only implemented for finite `p`"
        with pytest.raises(NotImplementedError, match=message):
            stats.pmean([2], np.inf)


    @pytest.mark.parametrize(
        "p, weights, ref",
        [
            # import numpy as np
            # from mpmath import mp
            # mp.dps = 1000
            #
            # x = [50., 100., 300., 70., 25., 100., 2.]
            # w = [1., 2., 1., 1., 1., 2., 1.]
            #
            # x = np.asarray([mp.mpf(x_) for x_ in x])
            # w = np.asarray([mp.mpf(w_) for w_ in w])
            # p = mp.mpf(2e-16)
            #
            # res = np.average(x**p)**(1/p)
            # float(res) # 47.23984586445856
            #
            (2e-16, None, 47.23984586445856),
            #
            # res = np.average(x**p, weights=w)**(1/p)
            # float(res) # 55.80644618389137
            #
            (2e-16, [1., 2., 1., 1., 1., 2., 1.], 55.80644618389137),
        ],
    )
    def test_small_p_gh23407(self, xp, p, weights, ref):
        # gh-23407 reported that pmean had
        # poor numerical behavior for very small nonzero p
        x = xp.asarray(