Port GeneralizedKSineSkewedWrappedCauchyDistribution (GSSC) from libDirectional

pyrecest/distributions/__init__.py

@@ -82,6 +82,7 @@
 from .circle.sine_skewed_distributions import (
     AbstractSineSkewedDistribution,
     GeneralizedKSineSkewedVonMisesDistribution,
+    GeneralizedKSineSkewedWrappedCauchyDistribution,
     GSSVMDistribution,
     SineSkewedVonMisesDistribution,
     SineSkewedWrappedCauchyDistribution,
@@ -275,6 +276,7 @@
     "CustomCircularDistribution",
     "AbstractSineSkewedDistribution",
     "GeneralizedKSineSkewedVonMisesDistribution",
+    "GeneralizedKSineSkewedWrappedCauchyDistribution",
     "GSSVMDistribution",
     "SineSkewedVonMisesDistribution",
     "SineSkewedWrappedCauchyDistribution",

pyrecest/distributions/circle/sine_skewed_distributions.py

@@ -1,7 +1,7 @@
 from abc import abstractmethod
 
 # pylint: disable=no-name-in-module,no-member
-from pyrecest.backend import mod, ndim, pi, sin
+from pyrecest.backend import cos, cosh, exp, mod, ndim, pi, sin, sinh
 from scipy.special import ive  # pylint: disable=no-name-in-module
 
 from .abstract_circular_distribution import AbstractCircularDistribution
@@ -180,3 +180,67 @@ def gamma(self):
 
     def base_pdf(self, xs):
         return self.wrapped_cauchy.pdf(xs)
+
+
+class GeneralizedKSineSkewedWrappedCauchyDistribution(AbstractCircularDistribution):
+    """
+    Generalized K Sine-Skewed Wrapped Cauchy (GSSC) distribution.
+    See Bekker, A., Nakhaei Rad, N., Arashi, M., Ley, C. (2020). Generalized Skew-Symmetric Circular and
+    Toroidal Distributions, Florence Nightingale Directional Statistics volume, Springer.
+    Parameters:
+    - mu (float): Mean direction of the distribution.
+    - gamma (float): Concentration parameter of the wrapped Cauchy distribution (positive).
+    - lambda_ (float): Skewness parameter, must be between -1 and 1 inclusive.
+    - k (int): Sine multiplier, currently supports only k=1.
+    - m (int): Power of the sine term, must be a positive integer.
+    """
+
+    # pylint: disable=too-many-positional-arguments
+    def __init__(self, mu, gamma, lambda_, k, m):
+        AbstractCircularDistribution.__init__(self)
+        self.mu = mod(mu, 2 * pi)
+        self.gamma = gamma
+        self.lambda_ = lambda_
+        self.k = k
+        self.m = m
+
+        self.validate_parameters()
+
+    def validate_parameters(self):
+        assert self.gamma > 0
+        assert -1.0 <= self.lambda_ and self.lambda_ <= 1.0
+        assert isinstance(self.m, int) and self.m >= 1
+
+    def pdf(self, xs):
+        assert self.k == 1, "Currently, only k=1 is supported"
+        # Use the WC pdf formula directly to ensure correct centering at mu
+        wc_pdf_vals = (
+            1 / (2 * pi) * sinh(self.gamma) / (cosh(self.gamma) - cos(xs - self.mu))
+        )
+        skew_factor = (1 + self.lambda_ * sin(self.k * (xs - self.mu))) ** self.m
+        # For the wrapped Cauchy: E[cos(n*(x-mu))] = exp(-n*k*gamma)
+        r2 = exp(-2 * self.k * self.gamma)
+        r4 = exp(-4 * self.k * self.gamma)
+        if self.m == 1:
+            norm_const = 1
+        elif self.m == 2:
+            norm_const = 1 / (1 + self.lambda_**2 / 2 * (1 - r2))
+        elif self.m == 3:
+            norm_const = 1 / (1 + 3 * self.lambda_**2 / 2 * (1 - r2))
+        elif self.m == 4:
+            norm_const = 1 / (
+                1
+                + self.lambda_**4 / 8 * (3 - 4 * r2 + r4)
+                + 3 * self.lambda_**2 * (1 - r2)
+            )
+        else:
+            raise NotImplementedError("m > 4 not implemented")
+
+        return norm_const * wc_pdf_vals * skew_factor
+
+    def shift(self, shift_by):
+        if ndim(shift_by) != 0:
+            raise ValueError("angle must be a scalar")
+        return GeneralizedKSineSkewedWrappedCauchyDistribution(
+            self.mu + shift_by, self.gamma, self.lambda_, self.k, self.m
+        )

pyrecest/tests/distributions/test_sine_skewed_distributions.py

@@ -1,8 +1,11 @@
 import unittest
 
+import numpy.testing as npt
+import pyrecest.backend
 from pyrecest.backend import array, pi
 from pyrecest.distributions.circle.sine_skewed_distributions import (
     GeneralizedKSineSkewedVonMisesDistribution,
+    GeneralizedKSineSkewedWrappedCauchyDistribution,
     SineSkewedWrappedCauchyDistribution,
     SineSkewedWrappedNormalDistribution,
 )
@@ -109,5 +112,99 @@ def test_sine_skewed_effect():
     assert skewed_dist.pdf(mu - 0.1) > normal_dist.pdf(mu)
 
 
+
+class TestGeneralizedKSineSkewedWrappedCauchyDistribution(unittest.TestCase):
+    def test_initialization(self):
+        """Test initialization with valid and invalid parameters."""
+        GeneralizedKSineSkewedWrappedCauchyDistribution(
+            mu=pi, gamma=0.5, lambda_=0.5, k=1, m=1
+        )
+
+        with self.assertRaises(AssertionError):
+            GeneralizedKSineSkewedWrappedCauchyDistribution(
+                mu=pi, gamma=0.5, lambda_=1.5, k=1, m=1
+            )
+
+        with self.assertRaises(AssertionError):
+            GeneralizedKSineSkewedWrappedCauchyDistribution(
+                mu=pi, gamma=0.5, lambda_=0.5, k=1, m=0
+            )
+
+        with self.assertRaises(AssertionError):
+            GeneralizedKSineSkewedWrappedCauchyDistribution(
+                mu=pi, gamma=-0.1, lambda_=0.5, k=1, m=1
+            )
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),  # pylint: disable=no-member
+        reason="Not supported on this backend",
+    )
+    def test_pdf_m1_normalizes(self):
+        """m=1 GSSC should integrate to 1."""
+        dist = GeneralizedKSineSkewedWrappedCauchyDistribution(
+            mu=pi / 3, gamma=0.5, lambda_=0.4, k=1, m=1
+        )
+        integral = dist.integrate_numerically()
+
+        npt.assert_allclose(integral, 1.0, atol=1e-4)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),  # pylint: disable=no-member
+        reason="Not supported on this backend",
+    )
+    def test_pdf_m2_normalizes(self):
+        """Test that m=2 PDF integrates to approximately 1."""
+        dist = GeneralizedKSineSkewedWrappedCauchyDistribution(
+            mu=pi / 4, gamma=0.3, lambda_=0.5, k=1, m=2
+        )
+        integral = dist.integrate_numerically()
+
+        npt.assert_allclose(integral, 1.0, atol=1e-4)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),  # pylint: disable=no-member
+        reason="Not supported on this backend",
+    )
+    def test_pdf_m3_normalizes(self):
+        """Test that m=3 PDF integrates to approximately 1."""
+        dist = GeneralizedKSineSkewedWrappedCauchyDistribution(
+            mu=0.0, gamma=0.5, lambda_=0.3, k=1, m=3
+        )
+        integral = dist.integrate_numerically()
+
+        npt.assert_allclose(integral, 1.0, atol=1e-4)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),  # pylint: disable=no-member
+        reason="Not supported on this backend",
+    )
+    def test_pdf_m4_normalizes(self):
+        """Test that m=4 PDF integrates to approximately 1."""
+        dist = GeneralizedKSineSkewedWrappedCauchyDistribution(
+            mu=pi, gamma=0.4, lambda_=0.6, k=1, m=4
+        )
+        integral = dist.integrate_numerically()
+
+        npt.assert_allclose(integral, 1.0, atol=1e-4)
+
+    def test_pdf_nonnegative(self):
+        """Test that PDF values are non-negative."""
+        for m in (1, 2, 3, 4):
+            dist = GeneralizedKSineSkewedWrappedCauchyDistribution(
+                mu=pi / 2, gamma=0.3, lambda_=0.5, k=1, m=m
+            )
+            xs = array([0.0, pi / 4, pi / 2, pi, 3 * pi / 2, 2 * pi - 0.01])
+            assert all(dist.pdf(xs) >= 0), f"Negative PDF value for m={m}"
+
+    def test_shift(self):
+        """Test the shift method modifies mu correctly."""
+        dist = GeneralizedKSineSkewedWrappedCauchyDistribution(
+            mu=0, gamma=0.5, lambda_=0.4, k=1, m=1
+        )
+        new_dist = dist.shift(array(pi / 2))
+
+        npt.assert_allclose(float(new_dist.mu), pi / 2, atol=1e-10)
+
+
 if __name__ == "__main__":
     unittest.main()