exact L-moments for gumbel_l

lmo/distributions/_genlambda.py

@@ -17,7 +17,7 @@
 import lmo.typing as lmt
 from lmo.special import harmonic
 from lmo.theoretical import entropy_from_qdf, l_moment_from_ppf
-from ._lm import get_lm_func
+from ._lm import lm_genlambda
 
 if sys.version_info >= (3, 13):
     from typing import override
@@ -129,7 +129,6 @@ def _genlambda_cdf0(  # noqa: C901
     [float],
     excluded={"ptol", "xtol", "maxiter"},
 )
-_genlambda_lm: Final = get_lm_func("genlambda")
 
 
 class genlambda_gen(rv_continuous):
@@ -187,7 +186,7 @@ def _stats(
         a, c = 1 + f, 1 - f
         b1, d1 = 1 + b, 1 + d
 
-        m1: float = 0 if b == d and f == 0 else _genlambda_lm(1, 0, 0, b, d, f).item()
+        m1 = 0 if b == d and f == 0 else float(lm_genlambda(1, 0, 0, b, d, f).item())
 
         if b <= -1 / 2 or d <= -1 / 2:
             return m1, math.nan, math.nan, math.nan
@@ -250,4 +249,4 @@ def _l_moment(
             )
             return np.asarray(lmbda_r)
 
-        return _genlambda_lm(r, s, t, b, d, f)
+        return lm_genlambda(r, s, t, b, d, f)

lmo/distributions/_lm.py

@@ -11,21 +11,31 @@
 import sys
 from collections.abc import Callable
 from math import gamma, log
-from typing import Any, Concatenate, Final, Literal, ParamSpec, TypeAlias, cast
-
-import numpy as np
-import optype.numpy as onp
-import scipy.special as sps
-
-import lmo.typing as lmt
-from lmo.special import harmonic
-from lmo.theoretical import l_moment_from_ppf
+from typing import (
+    TYPE_CHECKING,
+    Concatenate,
+    Final,
+    Literal,
+    ParamSpec,
+    Protocol,
+    TypeAlias,
+    TypeVar,
+    overload,
+)
 
 if sys.version_info >= (3, 13):
     from typing import TypeIs
 else:
     from typing_extensions import TypeIs
 
+import numpy as np
+import scipy.special as sps
+
+if TYPE_CHECKING:
+    import optype.numpy as onp
+
+from lmo.special import harmonic
+from lmo.theoretical import l_moment_from_ppf
 
 __all__ = [
     "lm_expon",
@@ -41,10 +51,11 @@
 
 DistributionName: TypeAlias = Literal[
     # 0 params
-    "uniform",
-    "logistic",
     "expon",
+    "gumbel_l",
     "gumbel_r",
+    "logistic",
+    "uniform",
     # 1 param
     "genextreme",
     "genpareto",
@@ -81,27 +92,47 @@
     "weibull_max",
 }
 
-_ArrF8: TypeAlias = onp.ArrayND[np.float64]
 
 _Tss = ParamSpec("_Tss")
+
 # (r, s, t, *params) -> float
 _LmFunc: TypeAlias = Callable[Concatenate[int, float, float, _Tss], float | np.float64]
 
-_LN2: Final = np.log(2)
-_LN3: Final = np.log(3)
-_LN5: Final = np.log(5)
+_F1_co = TypeVar("_F1_co", bound=_LmFunc[...], covariant=True)
 
 
-_LmVFunc: TypeAlias = lmt.Callable2[
-    Callable[
-        Concatenate[onp.ToIntND, float, float, _Tss],
-        onp.ArrayND[np.float64],
-    ],
-    Callable[
-        Concatenate[onp.ToInt, float, float, _Tss],
-        onp.Array[tuple[()], np.float64],
-    ],
-]
+class _LmVFuncWrapper(Protocol[_F1_co]):
+    @property
+    def pyfunc(self, /) -> _F1_co: ...
+
+    @overload
+    def __call__(
+        self: _LmVFuncWrapper[_LmFunc[_Tss]],
+        x: onp.ToInt,
+        s: float,
+        t: float,
+        /,
+        *args: _Tss.args,
+        **kwds: _Tss.kwargs,
+    ) -> onp.Array[tuple[()], np.float64]: ...
+    @overload
+    def __call__(
+        self: _LmVFuncWrapper[_LmFunc[_Tss]],
+        x: onp.ToIntND,
+        s: float,
+        t: float,
+        /,
+        *args: _Tss.args,
+        **kwds: _Tss.kwargs,
+    ) -> onp.ArrayND[np.float64]: ...
+
+
+_LmVFunc: TypeAlias = _LmVFuncWrapper[_LmFunc[_Tss]]
+
+
+_LN2: Final = log(2)
+_LN3: Final = log(3)
+_LN5: Final = log(5)
 
 
 def register_lm_func(
@@ -275,7 +306,14 @@ def lm_gumbel_r(r: int, s: float, t: float, /) -> np.float64 | float:
             return (_LN2 * -107 + _LN3 * 6 + _LN5 * 42) * 5 / 4
 
         case _:
-            return cast(Any, lm_genextreme).pyfunc(r, s, t, 0)
+            return lm_genextreme.pyfunc(r, s, t, 0)
+
+
+@register_lm_func("gumbel_l")
+def lm_gumbel_l(r: int, s: float, t: float, /) -> np.float64 | float:
+    """The reflected version of `lm_gumbel_r`."""
+    sign = -1 if r % 2 else 1
+    return sign * lm_gumbel_r.pyfunc(r, t, s)
 
 
 @register_lm_func("genextreme")
@@ -351,9 +389,9 @@ def lm_genpareto(r: int, s: float, t: float, /, a: float) -> np.float64 | float:
     if r == 0:
         return 1
     if a == 0:
-        return cast(Any, lm_expon).pyfunc(r, s, t)
+        return lm_expon.pyfunc(r, s, t)
     if a == 1:
-        return cast(Any, lm_uniform).pyfunc(r, s, t)
+        return lm_uniform.pyfunc(r, s, t)
 
     if not isinstance(t, int):
         msg = "fractional trimming"
@@ -399,9 +437,9 @@ def lm_kumaraswamy(
     return (
         np.sum(
             (-1) ** (k - t - 1)
-            * cast(_ArrF8, sps.comb(r + k - 2, r + t - 1))
-            * cast(_ArrF8, sps.comb(r + s + t, k))
-            * cast(_ArrF8, sps.beta(1 / a, 1 + k * b))
+            * sps.comb(r + k - 2, r + t - 1)
+            * sps.comb(r + s + t, k)
+            * sps.beta(1 / a, 1 + k * b)
             / a,
         )
         / r

tests/test_distributions.py

@@ -192,10 +192,11 @@ def test_genlambda(b: float, d: float, f: float):
 @pytest.mark.parametrize(
     "rv",
     [
-        distributions.uniform(),
-        distributions.logistic(),
-        distributions.expon(),
-        distributions.gumbel_r(),
+        distributions.uniform,
+        distributions.logistic,
+        distributions.expon,
+        distributions.gumbel_r,
+        distributions.gumbel_l,
         distributions.genextreme(-0.1),
         distributions.genpareto(0.1),
         kumaraswamy(2, 5),