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

Code before:
        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)


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)

    def test_levy_sf(self):
        # Large values, far into the tail of the distribution.
        x = np.array([1e15, 1e25, 1e35, 1e50])
        # Expected values were calculated with mpmath.
        expected = np.array([2.5231325220201597e-08,
                             2.52313252202016e-13,
                             2.52313252202016e-18,
                             7.978845608028653e-26])
        y = stats.levy.sf(x)
        assert_allclose(y, expected, rtol=1e-14)

    # The expected values for levy.isf(p) were calculated with mpmath.
    # For loc=0 and scale=1, the inverse SF can be computed with
    #
    #     import mpmath

Review comment (bug):
When you update `xsf` and remove the `xfail`, please choose either `assert_array_less` or `assert` for both. You can use `assert_array_less` here, if desired, by negating the first argument.

Code after:
        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])

        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)