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

Code before:
    def test_moment_gh18071(self, case):
        # gh-18071 reported an IntegrationWarning emitted by johnsonsu.stats
        # Check that the warning is no longer emitted and that the values
        # are accurate compared against results from Mathematica.
        # Reference values from Mathematica, e.g.
        # Mean[JohnsonDistribution["SU",-0.01, 1.1, 0.02, 0.0001]]
        res = stats.johnsonsu.stats(*case[:4], moments='mvsk')
        assert_allclose(res, case[4:], rtol=1e-14)


class TestTruncPareto:
    def test_pdf(self):
        # PDF is that of the truncated pareto distribution
        b, c = 1.8, 5.3
        x = np.linspace(1.8, 5.3)
        res = stats.truncpareto(b, c).pdf(x)
        ref = stats.pareto(b).pdf(x) / stats.pareto(b).cdf(c)
        assert_allclose(res, ref)

    def test_pdf_negative(self):
        # PDF is that of the truncated pareto distribution
        b, c = -1.8, 5.3
        x = np.linspace(1.8, 5.3)
        res = stats.truncpareto(b, c).pdf(x)
        ref = stats.pareto(b).pdf(x) / stats.pareto(b).cdf(c)
        assert_allclose(res, ref)

    @pytest.mark.parametrize('fix_loc', [True, False])
    @pytest.mark.parametrize('fix_scale', [True, False])
    @pytest.mark.parametrize('fix_b', [True, False])
    @pytest.mark.parametrize('fix_c', [True, False])
    def test_fit(self, fix_loc, fix_scale, fix_b, fix_c):

        rng = np.random.default_rng(6747363148258237171)
        b, c, loc, scale = 1.8, 5.3, 1, 2.5
        dist = stats.truncpareto(b, c, loc=loc, scale=scale)
        data = dist.rvs(size=500, random_state=rng)

        kwds = {}
        if fix_loc:
            kwds['floc'] = loc
        if fix_scale:
            kwds['fscale'] = scale
        if fix_b:
            kwds['f0'] = b
        if fix_c:
            kwds['f1'] = c

        if fix_loc and fix_scale and fix_b and fix_c:
            message = "All parameters fixed. There is nothing to optimize."
            with pytest.raises(RuntimeError, match=message):

Review comment (suggestion):
How about:
```suggestion
        # truncpareto is equivalent to more general powerlaw from gh-23648
        # exponent of truncpareto is negative in this case
        a, xmin, xmax = 4, 3, 5
        x = np.linspace(xmin, xmax)

        # compute reference using PDF from gh-23648
        C = a / (xmax ** a - xmin ** a)
        ref = C * x ** (a - 1)

        # compute using `truncpareto` with negative exponent
        b = -a
        c = xmax / xmin
        scale = xmin
        loc = 0
        res = stats.truncpareto(b, c, loc, scale).pdf(x)

        assert_allclose(res, ref)
```
This is what I meant when in https://github.com/scipy/scipy/issues/23648#issuecomment-3313762764 I wrote "where the reference value is obtained using your formulation of the powerlaw PDF."

Code after:
        """Tests fit for cases where floc is set.

        `rel_breitwigner` has special handling for these cases.
        """
        seed = 6936804688480013683
        rng = np.random.default_rng(seed)
        data = stats.rel_breitwigner.rvs(
            rho, scale=gamma, size=1000, random_state=rng
        )
        fit = stats.rel_breitwigner.fit(data, floc=0)
        assert_allclose((fit[0], fit[2]), (rho, gamma), rtol=2e-1)
        assert fit[1] == 0
        # Check again with fscale set.
        fit = stats.rel_breitwigner.fit(data, floc=0, fscale=gamma)
        assert_allclose(fit[0], rho, rtol=1e-2)
        assert (fit[1], fit[2]) == (0, gamma)


class TestJohnsonSU:
    @pytest.mark.parametrize("case", [  # a, b, loc, scale, m1, m2, g1, g2
            (-0.01, 1.1, 0.02, 0.0001, 0.02000137427557091,
             2.1112742956578063e-08, 0.05989781342460999, 20.36324408592951-3),
            (2.554395574161155, 2.2482281679651965, 0, 1, -1.54215386737391,
             0.7629882028469993, -1.256656139406788, 6.303058419339775-3)])
    def test_moment_gh18071(self, case):
        # gh-18071 reported an IntegrationWarning emitted by johnsonsu.stats
        # Check that the warning is no longer emitted and that the values
        # are accurate compared against results from Mathematica.
        # Reference values from Mathematica, e.g.
        # Mean[JohnsonDistribution["SU",-0.01, 1.1, 0.02, 0.0001]]
        res = stats.johnsonsu.stats(*case[:4], moments='mvsk')
        assert_allclose(res, case[4:], rtol=1e-14)


class TestTruncPareto:
    def test_pdf(self):
        # PDF is that of the truncated pareto distribution
        b, c = 1.8, 5.3
        x = np.linspace(1.8, 5.3)
        res = stats.truncpareto(b, c).pdf(x)
        ref = stats.pareto(b).pdf(x) / stats.pareto(b).cdf(c)
        assert_allclose(res, ref)

    def test_pdf_negative(self):
        # truncpareto is equivalent to more general powerlaw from gh-23648
        # exponent of truncpareto is negative in this case
        a, xmin, xmax = 4, 3, 5
        x = np.linspace(xmin, xmax)

        # compute reference using PDF from gh-23648
        C = a / (xmax ** a - xmin ** a)