diff --git a/src/pyrecest/distributions/circle/wrapped_laplace_distribution.py b/src/pyrecest/distributions/circle/wrapped_laplace_distribution.py index 34e63e79f..c5aceac4f 100644 --- a/src/pyrecest/distributions/circle/wrapped_laplace_distribution.py +++ b/src/pyrecest/distributions/circle/wrapped_laplace_distribution.py @@ -52,6 +52,21 @@ def _wrapped_exponential_density(rate, distance): return normalization * exp(-rate * distance) +def _mix_skew_components(positive_component, negative_component, kappa): + """Combine skew components without squaring a large ``kappa``.""" + if bool(all(kappa > 1.0)): + inverse_kappa = 1.0 / kappa + positive_weight = inverse_kappa * inverse_kappa + return (positive_weight * positive_component + negative_component) / ( + positive_weight + 1.0 + ) + + negative_weight = kappa * kappa + return (positive_component + negative_weight * negative_component) / ( + 1.0 + negative_weight + ) + + class WrappedLaplaceDistribution(AbstractCircularDistribution): """Wrapped Laplace distribution on the circle. @@ -86,9 +101,12 @@ def pdf(self, xs): xs = mod(xs, 2.0 * pi) positive_rate = self.lambda_ * self.kappa negative_rate = self.lambda_ / self.kappa - mixture_normalization = 1.0 + self.kappa**2 - p = ( - _wrapped_exponential_density(positive_rate, xs) - + self.kappa**2 * _wrapped_exponential_density(negative_rate, 2.0 * pi - xs) - ) / mixture_normalization - return p + positive_component = _wrapped_exponential_density(positive_rate, xs) + negative_component = _wrapped_exponential_density( + negative_rate, 2.0 * pi - xs + ) + return _mix_skew_components( + positive_component, + negative_component, + self.kappa, + ) diff --git a/tests/distributions/test_wrapped_laplace_distribution.py b/tests/distributions/test_wrapped_laplace_distribution.py index 5693b58c6..352218d08 100644 --- a/tests/distributions/test_wrapped_laplace_distribution.py +++ b/tests/distributions/test_wrapped_laplace_distribution.py @@ -1,6 +1,7 @@ import math import unittest +import numpy as np import numpy.testing as npt # pylint: disable=no-name-in-module,no-member @@ -70,6 +71,23 @@ def test_pdf_avoids_overflow_for_concentrated_negative_tail(self): npt.assert_allclose(actual, expected, rtol=1e-6) + @unittest.skipUnless( + pyrecest.backend.__backend_name__ == "numpy", + reason="Strict NumPy floating-point handling is backend-specific", + ) + def test_pdf_avoids_overflow_for_extreme_skew(self): + distribution = WrappedLaplaceDistribution(1.0, 1.0e200) + + with np.errstate( + over="raise", + invalid="raise", + divide="raise", + under="ignore", + ): + actual = distribution.pdf(array([0.0, 1.0])) + + npt.assert_allclose(actual, 1.0 / (2.0 * math.pi), rtol=1e-12) + def test_pdf_accepts_scalar_and_list_inputs(self): npt.assert_allclose(self.wl.pdf(1.0), self.wl.pdf(array(1.0)), rtol=1e-6) npt.assert_allclose(