Stabilize wrapped Laplace PDFs for extreme skew

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,
+        )

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(