From ea8199fade84be23f3828472edf0945e706c8eaa Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Fri, 29 Dec 2023 22:53:52 +0100
Subject: [PATCH 1/8] remove unused method

---
 lmo/special.py | 76 --------------------------------------------------
 1 file changed, 76 deletions(-)

diff --git a/lmo/special.py b/lmo/special.py
index be11651b..d54e3758 100644
--- a/lmo/special.py
+++ b/lmo/special.py
@@ -4,7 +4,6 @@
     'fpow',
     'gamma2',
     'harmonic',
-    'eval_sh_jacobi',
     'norm_sh_jacobi',
     'fourier_jacobi',
 )
@@ -204,81 +203,6 @@ def harmonic(
 
     return _out[()] if np.isscalar(n) else _out
 
-@overload
-def eval_sh_jacobi(
-    n: int,
-    alpha: float,
-    beta: float,
-    x: AnyScalar,
-    out: None = ...,
-) -> float: ...
-
-@overload
-def eval_sh_jacobi(
-    n: int,
-    alpha: float,
-    beta: float,
-    x: AnyNDArray[np.generic],
-    out: npt.NDArray[np.float64] | None = ...,
-) -> npt.NDArray[np.float64]: ...
-
-@overload
-def eval_sh_jacobi(
-    n: int,
-    alpha: float,
-    beta: float,
-    x: npt.ArrayLike,
-    out: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]: ...
-
-@overload
-def eval_sh_jacobi(
-    n: int,
-    alpha: float,
-    beta: float,
-    x: npt.ArrayLike,
-    out: npt.NDArray[np.float64] | None = ...,
-) -> float | npt.NDArray[np.float64]: ...
-
-def eval_sh_jacobi(
-    n: int,
-    alpha: float,
-    beta: float,
-    x: npt.ArrayLike,
-    out: npt.NDArray[np.float64] | None = None,
-) -> float | npt.NDArray[np.float64]:
-    r"""
-    Evaluate the (correct) shifted Jacobi Polynomial
-    \( \shjacobi{n}{\alpha}{\beta}{u} \) on \( x \in [0, 1] \), i.e.
-    the Jacobi Polynomial with mapped argument as \( x \mapsto 2x - 1 \).
-
-    It is well known that the "shifted Legendre" polynomial is the Legendre
-    polynomial with \( 2x - 1 \) as mapped argument, which is correctly
-    implemented in
-    [`scipy.special.eval_sh_legendre`][scipy.special.eval_sh_legendre].
-    Any textbook on orthogonal
-    polynomials will tell you that the generalization of Legendre are the
-    Jacobi polynomials. Hence, the only valid interpretation of the
-    "shifted Jacobi polynomials", should be the analogue (homomorphism) of
-    the shifted Legendre polynomials.
-
-    However, [`scipy.special.eval_sh_jacobi`][scipy.special.eval_sh_jacobi]
-    are **not** the shifted Jacobi polynomials!.
-    Instead, the method evaluates the *generalized Gegenbauer polynomials*.
-    The Jacobi-, and Legendre polynomials are denoted
-    with a "P", which stands for "polynomial". In the `eval_sh_jacobi`
-    docstring, the notation \( G_n^{p,q} \) is used. Clearly, the "G" stands
-    for "Gegenbauer".
-    See [scipy/scipy#18988](https://github.com/scipy/scipy/issues/18988) for
-    the relevant issue.
-    """
-    if alpha == beta == 0:
-        return sc.eval_sh_legendre(n, x, out)  # type: ignore
-
-    y = 2 * np.asanyarray(x) - 1
-    return sc.eval_jacobi(n, alpha, beta, y, out)  # type: ignore
-
-
 @overload
 def norm_sh_jacobi(n: int, alpha: float, beta: float) -> float: ...
 

From 9e4561758600dbcb866efff4c9ba9da4c3dbee6a Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 01:00:32 +0100
Subject: [PATCH 2/8] add `_poly` Jacobi extrema-related methods

---
 lmo/_poly.py | 187 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 187 insertions(+)

diff --git a/lmo/_poly.py b/lmo/_poly.py
index d0f51a3f..409f5399 100644
--- a/lmo/_poly.py
+++ b/lmo/_poly.py
@@ -9,6 +9,9 @@
 
 __all__ = (
     'eval_sh_jacobi',
+    'peaks_jacobi',
+    'arg_extrema_jacobi',
+    'extrema_jacobi',
     'jacobi',
     'jacobi_series',
     'roots',
@@ -101,6 +104,190 @@ def eval_sh_jacobi(
     return scs.eval_jacobi(n, a, b, u)  # type: ignore
 
 
+def peaks_jacobi(n: int, a: float, b: float) -> npt.NDArray[np.float64]:
+    r"""
+    Finds the \( x \in [-1, 1] \) s.t.
+    \( /frac{\dd{\shjacobi{n}{a}{b}{x}}}{\dd{x}} = 0 \) of a Jacobi polynomial,
+    which includes the endpoints \( x \in \{-1, 1\} \). I.e. the locations of
+    the peaks.
+
+    The Jacobi polynomials with order \( n \) have \( n + 1 \) peaks.
+
+    Examples:
+        For \( n = 0 \) there is only one "peak", since
+        \( \jacobi{0}{a}{b}{x} = 1 \):
+
+        >>> peaks_jacobi(0, 0, 0)
+        array([0.])
+
+        The `0` is arbitrary; all \( x \in [0, 1] \) evaluate to the same
+        constant \( 1 \).
+
+        For \( n = 1 \), it is a positive linear function, so the peaks
+        are exactly the endpoints, and do not depend on \( a \) or \( b \):
+
+        >>> peaks_jacobi(1, 0, 0)
+        array([-1., 1.])
+        >>> peaks_jacobi(1, 3.14, -1 / 12)
+        array([-1., 1.])
+
+        For \( n > 1 \), the effects of the choices for \( a \) and \( b \)
+        become apparant, e.g. for \( n = 4 \):
+
+        >>> peaks_jacobi(4, 0, 0).round(5)
+        array([-1.     , -0.65465,  0.     ,  0.65465,  1.     ])
+        >>> peaks_jacobi(4, 0, 1).round(5)
+        array([-1.     , -0.50779,  0.1323 ,  0.70882,  1.     ])
+        >>> peaks_jacobi(4, 1, 0).round(5)
+        array([-1.     , -0.70882, -0.1323 ,  0.50779,  1.     ])
+        >>> peaks_jacobi(4, 1, 1).round(5)
+        array([-1.     , -0.57735,  0.     ,  0.57735,  1.     ])
+        >>> peaks_jacobi(4, 2.5, 2.5)
+        array([-1. , -0.5,  0. ,  0.5,  1. ])
+        >>> peaks_jacobi(4, 10, 10).round(5)
+        array([-1.     , -0.33333,  0.     ,  0.33333,  1.     ])
+    """
+    if n == 0:
+        # constant; any x is a "peak"; so take the "middle ground"
+        return np.array([0.])
+    if n == 1:
+        # linear; the peaks are only at the ends
+        return np.array([-1., 1.])
+
+    # otherwise, peaks are at the ends, and at the roots of the derivative
+    x = np.empty(n + 1)
+    x[0] = -1
+    x[1:-1] = scs.roots_jacobi(n - 1, a + 1, b + 1)[0]  # type: ignore
+    x[-1] = 1
+
+    return np.round(x, 15) + 0.0  # cleanup of numerical noise
+
+def arg_extrema_jacobi(n: int, a: float, b: float) -> tuple[float, float]:
+    r"""
+    Find the \( x \) of the minimum and maximum values of a Jacobi polynomial
+    on \( [-1, 1] \).
+
+    Note:
+        There can be multiple \( x \) that share the same extremum, but only
+        one of them is returned, which for \( n > 0 \) is the smallest (first)
+        one.
+
+    Examples:
+        For \( n = 1 \), the Jacobi polynomials are positive linear function
+        (i.e. a straight line), so the minimum and maximum are the left and
+        right endpoints of the domain.
+
+        >>> arg_extrema_jacobi(1, 0, 0)
+        (-1.0, 1.0)
+        >>> arg_extrema_jacobi(1, 3.14, -1 / 12)
+        (-1.0, 1.0)
+
+        The 2nd degree Jacobi polynomial is a positive parabola, with one
+        unique minimum, and maxima at \( -1 \) and/or \( 1 \).
+        When \( a == b \), the parabola is centered within the domain, and
+        has maxima at both \( x = -1 \) and \( x=1 \). For the sake of
+        simplicity, only one (the first) value is returned in such cases:
+
+        >>> arg_extrema_jacobi(2, 0, 0)
+        (0.0, -1.0)
+        >>> arg_extrema_jacobi(2, 42, 42)
+        (0.0, -1.0)
+
+        Conversely, when \( a \neq b \), the parabola is "shifted" so that
+        there is only one global maximum:
+
+        >>> arg_extrema_jacobi(2, 0, 1)
+        (0.2, -1.0)
+        >>> arg_extrema_jacobi(2, 1, 0)
+        (-0.2, 1.0)
+        >>> arg_extrema_jacobi(2, 10, 2)
+        (-0.5, 1.0)
+
+    """
+    x = peaks_jacobi(n, a, b)
+    p = eval_sh_jacobi(n, a, b, (x + 1) / 2)
+
+    return x[np.argmin(p)], x[np.argmax(p)]
+
+def extrema_jacobi(n: int, a: float, b: float) -> tuple[float, float]:
+    r"""
+    Find the global minimum and maximum values of a (shifted) Jacobi
+    polynomial on \( [-1, 1] \) (or equivalently \( [0, 1] \) if shifted).
+
+    Examples:
+        With \( n \), \( \jacobi{0}{a}{b}{x} = 1 \), so there is only one
+        "extremum":
+
+        >>> extrema_jacobi(0, 0, 0)
+        (1, 1)
+        >>> extrema_jacobi(0, 3.14, -1 / 12)
+        (1, 1)
+
+        With \( n = 1 \), the extrema are always at \( -1 \) and \( 1 \),
+        but their values depend on \( a \) and \( b \):
+
+        >>> extrema_jacobi(1, 0, 0)
+        (-1.0, 1.0)
+        >>> extrema_jacobi(1, 0, 1)
+        (-2.0, 1.0)
+        >>> extrema_jacobi(1, 1, 0)
+        (-1.0, 2.0)
+        >>> extrema_jacobi(1, 1, 1)
+        (-2.0, 2.0)
+
+        For \( n = 2 \) (a parabola), the relation between \( a, b \)
+        and the extrema isn't as obvious:
+
+        >>> extrema_jacobi(2, 0, 0)
+        (-0.5, 1.0)
+        >>> extrema_jacobi(2, 0, 4)
+        (-0.75, 15.0)
+        >>> extrema_jacobi(2, 4, 0)
+        (-0.75, 15.0)
+        >>> extrema_jacobi(2, 4, 4)
+        (-1.5, 15.0)
+
+        With \( n = 3 \), the extrema appear to behave very predictable:
+
+        >>> extrema_jacobi(3, 0, 0)
+        (-1.0, 1.0)
+        >>> extrema_jacobi(3, 0, 1)
+        (-4.0, 1.0)
+        >>> extrema_jacobi(3, 1, 0)
+        (-1.0, 4.0)
+        >>> extrema_jacobi(3, 1, 1)
+        (-4.0, 4.0)
+
+        However, if we keep \( a \) fixed at \( 0 \), and increase \( b \),
+        the plot-twist emerges:
+
+        >>> extrema_jacobi(3, 0, 2)
+        (-10.0, 1.0)
+        >>> extrema_jacobi(3, 0, 3)
+        (-20.0, 1.0)
+        >>> extrema_jacobi(3, 0, 4)
+        (-35.0, 1.13541...)
+        >>> extrema_jacobi(3, 0, 5)
+        (-56.0, 1.25241...)
+
+        Looking at the corresponding \( x \) can help to understand the
+        "movement" of the maximum.
+
+        >>> arg_extrema_jacobi(3, 0, 2)
+        (-1.0, 1.0)
+        >>> arg_extrema_jacobi(3, 0, 3)
+        (-1.0, 0.0)
+        >>> arg_extrema_jacobi(3, 0, 4)
+        (-1.0, 0.09449...)
+        >>> arg_extrema_jacobi(3, 0, 5)
+        (-1.0, 0.17287...)
+
+    """
+    x = peaks_jacobi(n, a, b)
+    p = eval_sh_jacobi(n, a, b, (x + 1) / 2)
+    return cast(float, np.min(p)), cast(float, np.max(p))
+
+
 def _jacobi_coefs(n: int, a: float, b: float) -> npt.NDArray[np.float64]:
     p_n: np.poly1d
     p_n = scs.jacobi(n, a, b)  # type: ignore [reportUnknownMemberType]

From 916e98bfd70ff9da31722b4027abcf5d009c55fc Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 02:07:07 +0100
Subject: [PATCH 3/8] remove unused methods

---
 lmo/_poly.py | 41 -----------------------------------------
 1 file changed, 41 deletions(-)

diff --git a/lmo/_poly.py b/lmo/_poly.py
index 409f5399..8d000a6f 100644
--- a/lmo/_poly.py
+++ b/lmo/_poly.py
@@ -15,9 +15,6 @@
     'jacobi',
     'jacobi_series',
     'roots',
-    'extrema',
-    'minima',
-    'maxima',
 )
 
 from typing import Any, TypeVar, cast, overload
@@ -369,41 +366,3 @@ def roots(
         return x[(x >= a) & (x <= b)]
 
     return x
-
-
-def integrate(p: PolySeries, /, a: float | None = None) -> PolySeries:
-    r"""Calculate the anti-derivative: $P(x) = \int_a^x p(u) \, du$."""
-    return p.integ(lbnd=p.domain[0] if a is None else a)
-
-
-def extrema(
-    p: PolySeries,
-    /,
-    outside: bool = False,
-) -> npt.NDArray[np.inexact[Any]]:
-    """Return the $x$ in the domain of $p$, where $p'(x) = 0$."""
-    return roots(p.deriv(), outside=outside)
-
-
-def minima(
-    p: PolySeries,
-    /,
-    outside: bool = False,
-) -> npt.NDArray[np.inexact[Any]]:
-    """
-    Return the $x$ in the domain of $p$, where $p'(x) = 0$ and $p''(x) > 0$.
-    """  # noqa: D200
-    x = extrema(p, outside=outside)
-    return x[p.deriv(2)(x) > 0] if len(x) else x
-
-
-def maxima(
-    p: PolySeries,
-    /,
-    outside: bool = False,
-) -> npt.NDArray[np.inexact[Any]]:
-    """
-    Return the $x$ in the domain of $p$, where $p'(x) = 0$ and $p''(x) < 0$.
-    """  # noqa: D200
-    x = extrema(p, outside=outside)
-    return x[p.deriv(2)(x) < 0] if len(x) else x

From b33863b706484180cd37995df69480532d84439e Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 05:13:11 +0100
Subject: [PATCH 4/8] unused import

---
 lmo/_poly.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lmo/_poly.py b/lmo/_poly.py
index 8d000a6f..8d862363 100644
--- a/lmo/_poly.py
+++ b/lmo/_poly.py
@@ -17,7 +17,7 @@
     'roots',
 )
 
-from typing import Any, TypeVar, cast, overload
+from typing import TypeVar, cast, overload
 
 import numpy as np
 import numpy.polynomial as npp

From cf621ada1a7b528031144123c45f13ff8afc1d0d Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 05:18:50 +0100
Subject: [PATCH 5/8] significantly improved Hosking's L-ratio bounds

---
 lmo/diagnostic.py | 220 +++++++++++++++++++++++++---------------------
 1 file changed, 118 insertions(+), 102 deletions(-)

diff --git a/lmo/diagnostic.py b/lmo/diagnostic.py
index b260786b..efd681d5 100644
--- a/lmo/diagnostic.py
+++ b/lmo/diagnostic.py
@@ -41,7 +41,10 @@
     rv_frozen,
 )
 
+from .special import fpow
+
 from ._lm import l_ratio
+from ._poly import extrema_jacobi
 from ._utils import clean_orders, clean_trim
 from .typing import AnyInt, AnyTrim, IntVector
 
@@ -52,6 +55,8 @@
 
 AnyRV: TypeAlias = rv_continuous | rv_discrete
 
+_ArrF8: TypeAlias = npt.NDArray[np.float64]
+
 
 class HypothesisTestResult(NamedTuple):
     r"""
@@ -66,8 +71,8 @@ class HypothesisTestResult(NamedTuple):
             hypothesis, $H_0$.
     """
 
-    statistic: float | npt.NDArray[np.float64]
-    pvalue: float | npt.NDArray[np.float64]
+    statistic: float | _ArrF8
+    pvalue: float | _ArrF8
 
     @property
     def is_valid(self) -> bool | npt.NDArray[np.bool_]:
@@ -165,22 +170,14 @@ def normaltest(
     return HypothesisTestResult(k2, p_value)
 
 
-def _gof_stat_single(
-    l_obs: npt.NDArray[np.float64],
-    l_exp: npt.NDArray[np.float64],
-    cov: npt.NDArray[np.float64],
-) -> float:
+def _gof_stat_single(l_obs: _ArrF8, l_exp: _ArrF8, cov: _ArrF8) -> float:
     err = l_obs - l_exp
     prec = np.linalg.inv(cov)  # precision matrix
     return cast(float, err.T @ prec @ err)
 
 
 _gof_stat = cast(
-    Callable[[
-        npt.NDArray[np.float64],
-        npt.NDArray[np.float64],
-        npt.NDArray[np.float64],
-    ], npt.NDArray[np.float64]],
+    Callable[[_ArrF8, _ArrF8, _ArrF8], _ArrF8],
     np.vectorize(
         _gof_stat_single,
         otypes=[float],
@@ -192,7 +189,7 @@ def _gof_stat_single(
 
 def l_moment_gof(
     rv_or_cdf: AnyRV | Callable[[float], float],
-    l_moments: npt.NDArray[np.float64],
+    l_moments: _ArrF8,
     n_obs: int,
     /,
     trim: AnyTrim = (0, 0),
@@ -282,13 +279,13 @@ def l_moment_gof(
         lambda_rr = l_moment_cov_from_cdf(cdf, n, trim, **kwargs)
 
     stat = n_obs * _gof_stat(l_r.T, lambda_r, lambda_rr).T[()]
-    pval = cast(float | npt.NDArray[np.float64], chdtrc(n, stat))
+    pval = cast(float | _ArrF8, chdtrc(n, stat))
     return HypothesisTestResult(stat, pval)
 
 
 def l_stats_gof(
     rv_or_cdf: AnyRV | Callable[[float], float],
-    l_stats: npt.NDArray[np.float64],
+    l_stats: _ArrF8,
     n_obs: int,
     /,
     trim: AnyTrim = (0, 0),
@@ -316,7 +313,7 @@ def l_stats_gof(
         tau_rr = l_stats_cov_from_cdf(cdf, n, trim, **kwargs)
 
     stat = n_obs * _gof_stat(t_r.T, tau_r, tau_rr).T[()]
-    pval = cast(float | npt.NDArray[np.float64], chdtrc(n, stat))
+    pval = cast(float | _ArrF8, chdtrc(n, stat))
     return HypothesisTestResult(stat, pval)
 
 
@@ -348,7 +345,7 @@ def _lm2_bounds_single(r: int, trim: tuple[float, float]) -> float:
             ) / (np.pi * 2 * r**2)
 
 _lm2_bounds = cast(
-    Callable[[IntVector, tuple[float, float]], npt.NDArray[np.float64]],
+    Callable[[IntVector, tuple[float, float]], _ArrF8],
     np.vectorize(
         _lm2_bounds_single,
         otypes=[float],
@@ -357,33 +354,32 @@ def _lm2_bounds_single(r: int, trim: tuple[float, float]) -> float:
     ),
 )
 
-
 @overload
 def l_moment_bounds(
-    r: IntVector,
+    r: AnyInt,
     /,
     trim: AnyTrim = ...,
     scale: float = ...,
-) -> npt.NDArray[np.float64]:
+) -> float:
     ...
 
 
 @overload
 def l_moment_bounds(
-    r: AnyInt,
+    r: IntVector,
     /,
     trim: AnyTrim = ...,
     scale: float = ...,
-) -> float:
+) -> _ArrF8:
     ...
 
 
 def l_moment_bounds(
-    r: IntVector | AnyInt,
+    r: AnyInt | IntVector,
     /,
     trim: AnyTrim = (0, 0),
     scale: float = 1.0,
-) -> float | npt.NDArray[np.float64]:
+) -> float | _ArrF8:
     r"""
     Returns the absolute upper bounds $L^{(s,t)}_r$ on L-moments
     $\lambda^{(s,t)}_r$, proportional to the scale $\sigma_X$ (standard
@@ -480,109 +476,129 @@ def l_moment_bounds(
     _trim = clean_trim(trim)
     return scale * np.sqrt(_lm2_bounds(_r, _trim))[()]
 
-
 @overload
-def l_ratio_bounds(
-    r: IntVector,
-    /,
-    trim: AnyTrim = ...,
-    *,
-    has_variance: bool = ...,
-) -> npt.NDArray[np.float64]:
-    ...
-
-
+def l_ratio_bounds(r: AnyInt, /, trim: AnyTrim = ...) -> float: ...
 @overload
-def l_ratio_bounds(
-    r: AnyInt,
-    /,
-    trim: AnyTrim = ...,
-    *,
-    has_variance: bool = ...,
-) -> float:
-    ...
-
+def l_ratio_bounds(r: IntVector, /, trim: AnyTrim = ...) -> _ArrF8: ...
 
 def l_ratio_bounds(
     r: IntVector | AnyInt,
     /,
     trim: AnyTrim = (0, 0),
     *,
-    has_variance: bool = True,
-) -> float | npt.NDArray[np.float64]:
+    legacy: bool = False,
+) -> tuple[float | _ArrF8, float | _ArrF8]:
     r"""
-    Returns the absolute upper bounds $T^{(s,t)}_r$ on L-moment ratio's
-    $\tau^{(s,t)}_r = \lambda^{(s,t)}_r / \lambda^{(s,t)}_r$, for $r \ge 2$.
-    So $\left| \tau^{(s,t)}_r(X) \right| \le T^{(s,t)}_r$, given that
-    $\mathrm{Var}[X] = \sigma^2_X$ exists.
-
-    If the variance of the distribution is not defined, e.g. in case of the
-    [Cauchy distribution](https://wikipedia.org/wiki/Cauchy_distribution),
-    this method will not work. In this case, the looser bounds from
-    Hosking (2007) can be used instead, by passing `has_variance=False`.
-
-    Examples:
-        Calculate the bounds for different degrees of trimming:
-
-        >>> l_ratio_bounds([1, 2, 3, 4])
-        array([       inf, 1.        , 0.77459667, 0.65465367])
-        >>> # the bounds for r=1,2 are the same for all trimmings; skip them
-        >>> l_ratio_bounds([3, 4], trim=(1, 1))
-        array([0.61475926, 0.4546206 ])
-        >>> l_ratio_bounds([3, 4], trim=(.5, 1.5))
-        array([0.65060005, 0.49736473])
-
-        In case of undefined variance, the bounds become a lot looser:
-
-        >>> l_ratio_bounds([3, 4], has_variance=False)
-        array([1., 1.])
-        >>> l_ratio_bounds([3, 4], trim=(1, 1), has_variance=False)
-        array([1.11111111, 1.25      ])
-        >>> l_ratio_bounds([3, 4], trim=(.5, 1.5), has_variance=False)
-        array([1.33333333, 1.71428571])
+    Unlike the standardized product-moments, the L-moment ratio's with
+    \( r \ge 2 \) are bounded above and below.
+
+    Specifically, Hosking derived in 2007 that
+
+    \[
+        | \tlratio{r}{s,t}| \le
+            \frac 2 r
+            \frac{\ffact{r + s + t}{r - 2}}{\ffact{r - 1 + s \wedge t}{r - 2}}
+            .
+    \]
+
+    But this derivation relies on unnecessarily loose Jacobi polynomial bounds.
+    If the actual min and max of the Jacobi polynomials are used instead,
+    the following (tighter) inequality is obtained:
+
+    \[
+        \frac{\dot{w}_r^{(s, t)}}{\dot{w}_2^{(s, t)}}
+        \min_{u \in [0, 1]} \left[ \shjacobi{r - 1}{t + 1}{s + 1}{u} \right]
+        \le
+        \tlratio{s, t}{r}
+        \le
+        \frac{\dot{w}_r^{(s, t)}}{\dot{w}_2^{(s, t)}}
+        \max_{0 \le u \le 1} \left[ \shjacobi{r - 1}{t + 1}{s + 1}{u} \right],
+    \]
+
+    where
+
+    \[
+        \dot{w}_r^{(s, t)} =
+            % \frac{r + s + t}{r (r - 1)}
+            % \frac{\B(r,\ r + s + t)}{\B(r + s,\ r + t)}.
+            \frac{\B(r - 1,\ r + s + t + 1)}{r \B(r + s,\ r + t)}.
+    \]
 
     Args:
-        r: Order of the L-moment ratio(s), as positive integer scalar or
-            array-like.
-        trim: Tuple of left- and right- trim-lengths, matching those of the
-            relevant L-moment ratio's.
-        has_variance:
-            Set to False if the distribution has undefined variance, in which
-            case the (looser) bounds from J.R.M. Hosking (2007) will be used.
+        r: Scalar or array-like with the L-moment ratio order(s).
+        trim: L-moment ratio trim-length(s).
 
     Returns:
-        Array or scalar with shape like $r$.
+        A 2-tuple with arrays or scalars, of the lower- and upper bounds.
 
     See Also:
         - [`l_ratio`][lmo.l_ratio]
         - [`l_ratio_se`][lmo.l_ratio_se]
+        - [`diagnostic.l_moment_bounds`][lmo.diagnostic.l_moment_bounds]
 
     References:
         - [J.R.M. Hosking (2007) - Some theory and practical uses of trimmed
         L-moments](https://doi.org/10.1016/j.jspi.2006.12.002)
 
+    Todo:
+        - Fix the use in `distributions.py`
+        - Document `legacy` kwarg
+        - Examples
+            - Compare Hosking's (legacy) bounds with these ones
+            - Compare with L-stats of distributions, e.g. GPD with k < -35/9,
+            like mentioned by Hosking (2007).
+        - Some notes on how to find the extrema (i.e. d/du P = 0)
+
     """
-    if has_variance:
-        return l_moment_bounds(r, trim) / l_moment_bounds(2, trim)
+    _r = clean_orders(r)
+    s, t = clean_trim(trim)
 
-    # otherwise, fall back to the (very) loose bounds from Hosking
-    _r = clean_orders(r, rmin=1)
-    _n = cast(int, np.max(_r)) + 1
+    t_min = np.empty(_r.shape)
+    t_max = np.empty(_r.shape)
+
+    _cache: dict[int, tuple[float, float]] = {}
+    for i, ri in np.ndenumerate(_r):
+        _ri = cast(int, ri)
+        if _ri in _cache:
+            t_min[i], t_max[i] = _cache[_ri]
+
+        if _ri == 1:
+            # L-loc / L-scale; unbounded
+            t_min[i], t_max[i] = -np.inf, np.inf
+        elif _ri in (0, 2):  # or s == t == 0:
+            t_min[i] = t_max[i] = 1
+        elif legacy:
+            t_absmax = (
+                2
+                * fpow(_ri + s + t, _ri - 2)
+                / fpow(_ri + min(s, t) - 1, _ri - 2)
+                / _ri
+            )
+            t_min[i] = -t_absmax
+            t_max[i] = t_absmax
+        else:
+            cr_c2 = 2 * (
+                np.exp(
+                    lgamma(_ri - 1)
+                    - np.log(_ri)
+                    + lgamma(s + 2)
+                    - lgamma(_ri + s)
+                    + lgamma(t + 2)
+                    - lgamma(_ri + t)
+                    + lgamma(_ri + s + t + 1)
+                    - lgamma(s + t + 3)  # noqa: COM812
+                )
+            )
 
-    s, t = clean_trim(trim)
+            p_min, p_max = extrema_jacobi(_ri - 2, t + 1, s + 1)
+            assert p_min < 0 < p_max, (p_min, p_max)
 
-    out = np.ones(_n)
-    if _n > 1:
-        # `L-loc / L-scale (= 1 / CV)` is unbounded
-        out[1] = np.inf
+            t_min[i] = cr_c2 * p_min
+            t_max[i] = cr_c2 * p_max
 
-    if _n > 3 and s and t:
-        # if not trimmed, then the bounds are simply 1
-        p, q = s + t, min(s, t)
-        for _k in range(3, _n):
-            out[_k] = out[_k - 1] * (1 + p / _k) / (1 + q / (_k - 1))
+        _cache[_ri] = t_min[i], t_max[i]
 
-    return out[_r]
+    return t_min[()], t_max[()]
 
 
 def rejection_point(
@@ -668,7 +684,7 @@ def rejection_point(
     def integrand(x: float) -> float:
         return max(abs(influence_fn(-x)), abs(influence_fn(x)))
 
-    def obj(r: npt.NDArray[np.float64]) -> float:
+    def obj(r: _ArrF8) -> float:
         return quad(integrand, r[0], np.inf)[0]  # type: ignore
 
     res = cast(
@@ -734,7 +750,7 @@ def error_sensitivity(
     if np.isinf(influence_fn(a)) or np.isinf(influence_fn(b)):
         return np.inf
 
-    def obj(xs: npt.NDArray[np.float64]) -> float:
+    def obj(xs: _ArrF8) -> float:
         return -abs(influence_fn(xs[0]))
 
     bounds = None if np.isneginf(a) and np.isposinf(b) else [(a, b)]
@@ -817,7 +833,7 @@ def shift_sensitivity(
 
     """
 
-    def obj(xs: npt.NDArray[np.float64]) -> float:
+    def obj(xs: _ArrF8) -> float:
         x, y = xs
         if y == x:
             return 0

From dcd081c832f2c7c4b57ff0c8b63a859e111c7c11 Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 21:03:14 +0100
Subject: [PATCH 6/8] fixed outdated `l_ratio_bounds` usages

---
 lmo/distributions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lmo/distributions.py b/lmo/distributions.py
index 34665332..02377e3c 100644
--- a/lmo/distributions.py
+++ b/lmo/distributions.py
@@ -848,7 +848,7 @@ def _check_lmoments(
 
     r = np.arange(1, n + 1)
     t_r = l_r[2:] / l_r[1]
-    t_r_max = l_ratio_bounds(r[2:], trim, has_variance=False)
+    t_r_max = l_ratio_bounds(r[2:], trim, legacy=True)[1]
     if np.any(rs0_oob := np.abs(t_r) > t_r_max):
         r_oob = np.argwhere(rs0_oob)[0] + 3
         t_oob = t_r[rs0_oob][0]

From 9749ae9b1e4b7fa31be7e316c66474ee1ae99d0a Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 21:03:43 +0100
Subject: [PATCH 7/8] doctests and finishing touches in `l_ratio_bounds`

---
 lmo/diagnostic.py | 75 ++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 61 insertions(+), 14 deletions(-)

diff --git a/lmo/diagnostic.py b/lmo/diagnostic.py
index efd681d5..83654c8c 100644
--- a/lmo/diagnostic.py
+++ b/lmo/diagnostic.py
@@ -41,11 +41,10 @@
     rv_frozen,
 )
 
-from .special import fpow
-
 from ._lm import l_ratio
 from ._poly import extrema_jacobi
 from ._utils import clean_orders, clean_trim
+from .special import fpow
 from .typing import AnyInt, AnyTrim, IntVector
 
 if TYPE_CHECKING:
@@ -477,9 +476,22 @@ def l_moment_bounds(
     return scale * np.sqrt(_lm2_bounds(_r, _trim))[()]
 
 @overload
-def l_ratio_bounds(r: AnyInt, /, trim: AnyTrim = ...) -> float: ...
+def l_ratio_bounds(
+    r: AnyInt,
+    /,
+    trim: AnyTrim = ...,
+    *,
+    legacy: bool = ...,
+) -> tuple[float, float]: ...
+
 @overload
-def l_ratio_bounds(r: IntVector, /, trim: AnyTrim = ...) -> _ArrF8: ...
+def l_ratio_bounds(
+    r: IntVector,
+    /,
+    trim: AnyTrim = ...,
+    *,
+    legacy: bool = ...,
+) -> tuple[_ArrF8, _ArrF8]: ...
 
 def l_ratio_bounds(
     r: IntVector | AnyInt,
@@ -524,9 +536,53 @@ def l_ratio_bounds(
             \frac{\B(r - 1,\ r + s + t + 1)}{r \B(r + s,\ r + t)}.
     \]
 
+    Examples:
+        Without trim, the lower- and upper-bounds of the L-skewness and
+        L-kurtosis are:
+
+        >>> l_ratio_bounds(3)
+        (-1.0, 1.0)
+        >>> l_ratio_bounds(4)
+        (-0.25, 1.0)
+
+        For the L-kurtosis, the "legacy" bounds by Hosking (2007) are clearly
+        looser:
+
+        >>> l_ratio_bounds(4, legacy=True)
+        (-1.0, 1.0)
+
+        For the symmetrically trimmed TL-moment ratio's:
+
+        >>> l_ratio_bounds(3, trim=3)
+        (-1.2, 1.2)
+        >>> l_ratio_bounds(4, trim=3)
+        (-0.15, 1.5)
+
+        Simiarly, those of the LL-ratio's are
+
+        >>> l_ratio_bounds(3, trim=(0, 3))
+        (-0.8, 2.0)
+        >>> l_ratio_bounds(4, trim=(0, 3))
+        (-0.2333..., 3.5)
+
+        The LH-skewness bounds are "flipped" w.r.t to the LL-skewness,
+        but they are the same for the L*-kurtosis:
+
+        >>> l_ratio_bounds(3, trim=(3, 0))
+        (-2.0, 0.8)
+        >>> l_ratio_bounds(4, trim=(3, 0))
+        (-0.2333..., 3.5)
+
+        The bounds of multiple L-ratio's can be calculated in one shot:
+        >>> l_ratio_bounds([3, 4, 5, 6], trim=(1, 2))
+        (array([-1.        , -0.194..., -1.12      , -0.149...]),
+         array([1.333..., 1.75      , 2.24      , 2.8       ]))
+
+
     Args:
         r: Scalar or array-like with the L-moment ratio order(s).
         trim: L-moment ratio trim-length(s).
+        legacy: If set to `True`, will use the (looser) by Hosking (2007).
 
     Returns:
         A 2-tuple with arrays or scalars, of the lower- and upper bounds.
@@ -540,15 +596,6 @@ def l_ratio_bounds(
         - [J.R.M. Hosking (2007) - Some theory and practical uses of trimmed
         L-moments](https://doi.org/10.1016/j.jspi.2006.12.002)
 
-    Todo:
-        - Fix the use in `distributions.py`
-        - Document `legacy` kwarg
-        - Examples
-            - Compare Hosking's (legacy) bounds with these ones
-            - Compare with L-stats of distributions, e.g. GPD with k < -35/9,
-            like mentioned by Hosking (2007).
-        - Some notes on how to find the extrema (i.e. d/du P = 0)
-
     """
     _r = clean_orders(r)
     s, t = clean_trim(trim)
@@ -598,7 +645,7 @@ def l_ratio_bounds(
 
         _cache[_ri] = t_min[i], t_max[i]
 
-    return t_min[()], t_max[()]
+    return t_min.round(12)[()], t_max.round(12)[()]
 
 
 def rejection_point(

From 413dc8c01f8c1c13ad8a671f89bc931fe7af44b0 Mon Sep 17 00:00:00 2001
From: jorenham <jhammudoglu@gmail.com>
Date: Sat, 30 Dec 2023 21:04:24 +0100
Subject: [PATCH 8/8] typo fixes

---
 lmo/_poly.py      | 2 +-
 lmo/diagnostic.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/lmo/_poly.py b/lmo/_poly.py
index 8d862363..0d59784f 100644
--- a/lmo/_poly.py
+++ b/lmo/_poly.py
@@ -129,7 +129,7 @@ def peaks_jacobi(n: int, a: float, b: float) -> npt.NDArray[np.float64]:
         array([-1., 1.])
 
         For \( n > 1 \), the effects of the choices for \( a \) and \( b \)
-        become apparant, e.g. for \( n = 4 \):
+        become apparent, e.g. for \( n = 4 \):
 
         >>> peaks_jacobi(4, 0, 0).round(5)
         array([-1.     , -0.65465,  0.     ,  0.65465,  1.     ])
diff --git a/lmo/diagnostic.py b/lmo/diagnostic.py
index 83654c8c..f43e6ce1 100644
--- a/lmo/diagnostic.py
+++ b/lmo/diagnostic.py
@@ -558,7 +558,7 @@ def l_ratio_bounds(
         >>> l_ratio_bounds(4, trim=3)
         (-0.15, 1.5)
 
-        Simiarly, those of the LL-ratio's are
+        Similarly, those of the LL-ratio's are
 
         >>> l_ratio_bounds(3, trim=(0, 3))
         (-0.8, 2.0)