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

Code before:
                return {'a': {'endpoints': (0, np.inf), 'inclusive': (False, False)},
                        'b': {'endpoints': ('a', np.inf), 'inclusive': (False, False)}}

            @property
            def support(self):
                return {'endpoints': ('a', 'b')}

            def pdf(self, x, a, b):
                return 1 / (x * (np.log(b) - np.log(a)))

            def sample(self, shape, *, a, b, rng=None):
                p = rng.uniform(size=shape)
                return np.exp(np.log(a) + p * (np.log(b) - np.log(a)))

            def moment(self, order, kind='raw', *, a, b):
                if order == 1 and kind == 'raw':
                    # quadrature is perfectly accurate here; add 1e-10 error so we
                    # can tell the difference between the two
                    return (b - a) / np.log(b/a) + 1e-10

            def lmoment(self, order, *, a, b):
                s = np.log(b/a)
                l1 = (b - a) / s
                l2 = ((s - 2) * b + (s + 2) * a) / s**2
                l3 = ((s**2 - 6*s + 12) * b - (s**2 + 6*s + 12) * a) / s ** 3
                l4 = ((s**3 - 12*s**2 + 60*s - 120) * b
                      + (s**3 + 12*s**2 + 60*s + 120) * a) / s ** 4
                ls = {1: l1, 2: l2, 3: l3, 4: l4}
                return ls.get(int(order), None)

        LogUniform = stats.make_distribution(MyLogUniform())

        X = LogUniform(a=1., b=np.e)
        Y = stats.exp(Uniform(a=0., b=1.))

        # pre-2.0 support is not needed for much longer, so let's just test with 2.0+
        if np.__version__ >= "2.0":
            assert str(X) == f"MyLogUniform(a=1.0, b={np.e})"
            assert repr(X) == f"MyLogUniform(a=np.float64(1.0), b=np.float64({np.e}))"

        x = X.sample(shape=10, rng=rng)
        p = X.cdf(x)

        assert_allclose(X.support(), Y.support())
        assert_allclose(X.entropy(), Y.entropy())
        assert_allclose(X.median(), Y.median())
        assert_allclose(X.logpdf(x), Y.logpdf(x))
        assert_allclose(X.pdf(x), Y.pdf(x))
        assert_allclose(X.logcdf(x), Y.logcdf(x))
        assert_allclose(X.cdf(x), Y.cdf(x))
        assert_allclose(X.logccdf(x), Y.logccdf(x))

Review comment (refactor):
To verify against definition from [here](https://en.wikipedia.org/wiki/L-moment):

<img width="966" height="270" alt="image" src="https://github.com/user-attachments/assets/959d0ba1-2af5-42b3-ad78-fc754c2b079a" />

```python3
from sympy import *
p, a, b, c, x = symbols('p, a, b, c, x', positive=True, real=True)

Q = solve( (log(x) - log(a))/(log(b)- log(a)) - p, x)[0]
l1_ref = simplify(integrate(Q, (p, 0, 1)))
l2_ref = simplify(integrate(Q*(2*p - 1), (p, 0, 1)))
l3_ref = simplify(integrate(Q*(6*p**2 - 6*p + 1), (p, 0, 1)))
l4_ref = simplify(integrate(Q*(20*p**3 -30*p**2 + 12*p - 1), (p, 0, 1)))

s = log(b/a)
l1 = (b - a)/s
simplify(l1 - l1_ref)

l2 = ((s - 2) * b + (s + 2) * a) / s**2
simplify(l2 - l2_ref)

l3 = ((s**2 - 6*s + 12) * b - (s**2 + 6*s + 12) * a) / s ** 3
simplify(l3 - l3_ref)

l4 = ((s**3 - 12*s**2 + 60*s - 120) * b + (s**3 + 12*s**2 + 60*s + 120) * a) / s ** 4
simplify(l4 - l4_ref)
```

Code after:
            @property
            def support(self):
                return {'endpoints': ('a', 'b')}

            def pdf(self, x, a, b):
                return 1 / (x * (np.log(b) - np.log(a)))

            def sample(self, shape, *, a, b, rng=None):
                p = rng.uniform(size=shape)
                return np.exp(np.log(a) + p * (np.log(b) - np.log(a)))

            def moment(self, order, kind='raw', *, a, b):
                if order == 1 and kind == 'raw':
                    # quadrature is perfectly accurate here; add 1e-10 error so we
                    # can tell the difference between the two
                    return (b - a) / np.log(b/a) + 1e-10

            def lmoment(self, order, *, a, b):
                s = np.log(b/a)
                l1 = (b - a) / s
                l2 = ((s - 2) * b + (s + 2) * a) / s**2
                l3 = ((s**2 - 6*s + 12) * b - (s**2 + 6*s + 12) * a) / s ** 3
                l4 = ((s**3 - 12*s**2 + 60*s - 120) * b
                      + (s**3 + 12*s**2 + 60*s + 120) * a) / s ** 4
                ls = {1: l1, 2: l2, 3: l3, 4: l4}
                return ls.get(int(order), None)

        LogUniform = stats.make_distribution(MyLogUniform())

        X = LogUniform(a=1., b=np.e)
        Y = stats.exp(Uniform(a=0., b=1.))

        assert str(X) == f"MyLogUniform(a=1.0, b={np.e})"
        assert repr(X) == f"MyLogUniform(a=np.float64(1.0), b=np.float64({np.e}))"

        x = X.sample(shape=10, rng=rng)
        p = X.cdf(x)

        assert_allclose(X.support(), Y.support())
        assert_allclose(X.entropy(), Y.entropy())
        assert_allclose(X.median(), Y.median())
        assert_allclose(X.logpdf(x), Y.logpdf(x))
        assert_allclose(X.pdf(x), Y.pdf(x))
        assert_allclose(X.logcdf(x), Y.logcdf(x))
        assert_allclose(X.cdf(x), Y.cdf(x))
        assert_allclose(X.logccdf(x), Y.logccdf(x))
        assert_allclose(X.ccdf(x), Y.ccdf(x))
        assert_allclose(X.icdf(p), Y.icdf(p))
        assert_allclose(X.iccdf(p), Y.iccdf(p))
        for kind in ['raw', 'central', 'standardized']:
            for order in range(5):