diff --git a/pyrecest/distributions/__init__.py b/pyrecest/distributions/__init__.py index 1f73ca5d1..e96b4549f 100644 --- a/pyrecest/distributions/__init__.py +++ b/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", diff --git a/pyrecest/distributions/circle/sine_skewed_distributions.py b/pyrecest/distributions/circle/sine_skewed_distributions.py index 1bb8eabd3..d61623ca4 100644 --- a/pyrecest/distributions/circle/sine_skewed_distributions.py +++ b/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 + ) diff --git a/pyrecest/tests/distributions/test_sine_skewed_distributions.py b/pyrecest/tests/distributions/test_sine_skewed_distributions.py index 7686e6483..ee96a6b8d 100644 --- a/pyrecest/tests/distributions/test_sine_skewed_distributions.py +++ b/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()