Repository: scipy/scipy
File: scipy/stats/_distribution_infrastructure.py (Python)

Code before:
        return lmoment

    def _lmoment_dispatch(self, order, *, methods, **params):
        lmoment = None

        if 'cache' in methods:
            lmoment = self._lmoment_cache.get(order, None)

        if lmoment is None and 'formula' in methods:
            lmoment = self._lmoment_formula(order, **params)

        if lmoment is None and 'general' in methods:
            lmoment = self._lmoment_general(order, **params)

        if lmoment is None and 'quadrature_icdf' in methods and (
                self._overrides('_icdf_formula') or self._overrides('_iccdf_formula')):
            lmoment = self._lmoment_integrate_icdf(order, **params)

        if lmoment is None and 'order_statistics' in methods:
            lmoment = self._lmoment_from_order_statistics(order, **params)

        if lmoment is None and 'quadrature_icdf' in methods:
            lmoment = self._lmoment_integrate_icdf(order, **params)

        if lmoment is not None and self.cache_policy != _NO_CACHE:
            self._lmoment_cache[order] = lmoment

        return lmoment

    def _lmoment_formula(self, order, **params):
        return None

    def _lmoment_general(self, order, **params):
        return self.mean() if order == 1 else None

    def _lmoment_from_order_statistics(self, order, **params):
        k = np.arange(order)
        k = xpx.atleast_nd(k, ndim=self._ndim + 1).T
        E = order_statistic(self, r=order-k, n=order).mean()
        bc = special.binom(order-1, k)
        return np.sum((-1)**k * bc * E, axis=0) / order

    def _lmoment_integrate_icdf(self, order, **params):
        def integrand(p, **params):
            x = self._icdf_dispatch(p, **params)
            P = special.eval_sh_legendre(order - 1, p)
            return x * P
        return self._quadrature(integrand, limits=(0., 1.), params=params)

    ### Convenience


Review comment (refactor):
It's also possible to use quadrature with the CDF (without the PDF).
From https://jorenham.github.io/Lmo/api/low_level/#lmo.theoretical.l_moment_from_cdf:

<img width="793" height="509" alt="image" src="https://github.com/user-attachments/assets/05f9f1d7-690d-4ec3-91e1-4d619b3b2ac5" />

This is for the generalized trimmed L-moments, so for the L-moments you set s=t=0, and the maths will simplify quite a bit. Note that the shifted Jacobi polynomial there isn't  the `eval_sh_jacobi`, but rather `eval_jacobi(r - 2, 1, 1, 2*u - 1)`.

[implementation](https://github.com/jorenham/Lmo/blob/38215eb525eb920e440c3e1363e223720f1d2ccb/lmo/theoretical/_f_to_lm.py#L216-L261)

Code after:
    def _lmoment_dispatch(self, order, *, methods, **params):
        lmoment = None

        if 'cache' in methods:
            lmoment = self._lmoment_cache.get(order, None)

        if lmoment is None and 'formula' in methods:
            lmoment = self._lmoment_formula(order, **params)

        if lmoment is None and 'general' in methods:
            lmoment = self._lmoment_general(order, **params)

        if lmoment is None and 'quadrature_icdf' in methods and (
                self._overrides('_icdf_formula') or self._overrides('_iccdf_formula')):
            lmoment = self._lmoment_integrate_icdf(order, **params)

        if lmoment is None and 'order_statistics' in methods:
            lmoment = self._lmoment_from_order_statistics(order, **params)

        if lmoment is None and 'quadrature_icdf' in methods:
            lmoment = self._lmoment_integrate_icdf(order, **params)

        if lmoment is not None and self.cache_policy != _NO_CACHE:
            self._lmoment_cache[order] = lmoment

        return lmoment

    def _lmoment_formula(self, order, **params):
        return None

    def _lmoment_general(self, order, **params):
        return self.mean() if order == 1 else None

    def _lmoment_from_order_statistics(self, order, **params):
        k = np.arange(order)
        k = xpx.atleast_nd(k, ndim=self._ndim + 1).T
        E = order_statistic(self, r=order-k, n=order).mean()
        bc = special.binom(order-1, k)
        return np.sum((-1)**k * bc * E, axis=0) / order

    def _lmoment_integrate_icdf(self, order, **params):
        def integrand(p, **params):
            x = self._icdf_dispatch(p, **params)
            P = special.eval_sh_legendre(order - 1, p)
            return x * P
        return self._quadrature(integrand, limits=(0., 1.), params=params)

    ### Convenience

    def plot(self, x='x', y=None, *, t=None, ax=None):
        r"""Plot a function of the distribution.