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

Code before:
        # See test_tukeylamdba_stats.py for more tests.

        mv = stats.tukeylambda.stats(0, moments='mvsk')
        # Known exact values:
        expected = [0, np.pi**2/3, 0, 1.2]
        assert_almost_equal(mv, expected, decimal=10)

        mv = stats.tukeylambda.stats(3.13, moments='mvsk')
        # 'expected' computed with mpmath.
        expected = [0, 0.0269220858861465102, 0, -0.898062386219224104]
        assert_almost_equal(mv, expected, decimal=10)

        mv = stats.tukeylambda.stats(0.14, moments='mvsk')
        # 'expected' computed with mpmath.
        expected = [0, 2.11029702221450250, 0, -0.02708377353223019456]
        assert_almost_equal(mv, expected, decimal=10)

    def test_large_argument_ticket_21370(self):
        # Check that survival function goes to 0 as x gets large
        x = [1e4, 1e6, 1e8, 1e10, 1e12, 1e24]
        lam = -0.5
        probs = stats.tukeylambda.sf(x, lam)
        # Check that the sf is 0 for large x
        assert_(probs[-1] == 0.0)
        # Check that the sf is always decreasing except when it is 0
        assert_(np.all(np.diff(probs)[probs[:-1] > 0.0]))
        
class TestLevy:

    def test_levy_cdf_ppf(self):
        # Test levy.cdf, including small arguments.
        x = np.array([1000, 1.0, 0.5, 0.1, 0.01, 0.001])

        # Expected values were calculated separately with mpmath.
        # E.g.
        # >>> mpmath.mp.dps = 100
        # >>> x = mpmath.mp.mpf('0.01')
        # >>> cdf = mpmath.erfc(mpmath.sqrt(1/(2*x)))
        expected = np.array([0.9747728793699604,
                             0.3173105078629141,
                             0.1572992070502851,
                             0.0015654022580025495,
                             1.523970604832105e-23,
                             1.795832784800726e-219])

        y = stats.levy.cdf(x)
        assert_allclose(y, expected, rtol=1e-10)

        # ppf(expected) should get us back to x.
        xx = stats.levy.ppf(expected)
        assert_allclose(xx, x, rtol=1e-13)

Review comment (suggestion):
```suggestion
        assert np.all(np.diff(probs)[probs[:-1] < 0.0])
```
Right? 

Also, there's the option to use https://numpy.org/doc/2.1/reference/generated/numpy.testing.assert_array_less.html.

Code after:
        # See test_tukeylamdba_stats.py for more tests.

        mv = stats.tukeylambda.stats(0, moments='mvsk')
        # Known exact values:
        expected = [0, np.pi**2/3, 0, 1.2]
        assert_almost_equal(mv, expected, decimal=10)

        mv = stats.tukeylambda.stats(3.13, moments='mvsk')
        # 'expected' computed with mpmath.
        expected = [0, 0.0269220858861465102, 0, -0.898062386219224104]
        assert_almost_equal(mv, expected, decimal=10)

        mv = stats.tukeylambda.stats(0.14, moments='mvsk')
        # 'expected' computed with mpmath.
        expected = [0, 2.11029702221450250, 0, -0.02708377353223019456]
        assert_almost_equal(mv, expected, decimal=10)

    @pytest.mark.xfail(reason="Precision in tail hasn't been fixed in xsf")
    def test_large_argument_ticket_21370(self):
        # Check that survival function goes to 0 as x gets large
        x = [1e4, 1e6, 1e8, 1e10, 1e12, 1e24]
        lam = -0.5
        sf = stats.tukeylambda.sf(x, lam)
        # Check that the sf is 0 for large x
        assert sf[-1] == 0.0
        # Check that the sf is always decreasing except when it is 0
        assert_array_less(np.diff(sf[sf > 0]), 0)

        cdf = stats.tukeylambda.cdf(x, lam)
        assert cdf[-1] == 1.0
        # Check that the cdf is always increasing except when it saturates 
        # to 1.0
        assert np.all(np.diff(cdf[cdf < 1]) > 0)
        
class TestLevy:

    def test_levy_cdf_ppf(self):
        # Test levy.cdf, including small arguments.
        x = np.array([1000, 1.0, 0.5, 0.1, 0.01, 0.001])

        # Expected values were calculated separately with mpmath.
        # E.g.
        # >>> mpmath.mp.dps = 100
        # >>> x = mpmath.mp.mpf('0.01')
        # >>> cdf = mpmath.erfc(mpmath.sqrt(1/(2*x)))
        expected = np.array([0.9747728793699604,
                             0.3173105078629141,
                             0.1572992070502851,
                             0.0015654022580025495,
                             1.523970604832105e-23,
                             1.795832784800726e-219])