Support negative-kappa Watson Bingham conversion

src/pyrecest/distributions/hypersphere_subset/watson_distribution.py

@@ -166,17 +166,25 @@ def ln_pdf(self, xs):
         return self.ln_norm_const + self.kappa * (xs @ self.mu) ** 2
 
     def to_bingham(self) -> BinghamDistribution:
-        if self.kappa < 0:
-            raise NotImplementedError(
-                "Conversion to Bingham is not implemented for kappa<0"
-            )
-
         M = tile(self.mu.reshape(-1, 1), (1, self.input_dim))
         E = diag(array(concatenate((array([0]), ones(self.input_dim - 1)))))
         M = M + E
         Q, _ = linalg.qr(M)
-        M = hstack([Q[:, 1:], Q[:, 0].reshape(-1, 1)])
-        Z = hstack((full((self.dim,), -self.kappa), array(0.0)))
+
+        if self.kappa >= 0:
+            # Bingham concentrations must be ascending with a final zero.  Put
+            # the Watson axis last and shift the exponent by -kappa, which does
+            # not change the normalized density.
+            M = hstack([Q[:, 1:], Q[:, 0].reshape(-1, 1)])
+            Z = hstack((full((self.dim,), -self.kappa), array(0.0)))
+        else:
+            # For negative kappa the Watson axis is the least likely direction.
+            # Keeping it as the first Bingham eigenvector gives the exponent
+            # kappa * (mu.T @ x)^2 directly, while the orthogonal complement
+            # receives zero concentration.
+            M = Q
+            Z = hstack((array([self.kappa]), zeros(self.dim)))
+
         return BinghamDistribution(Z, M)
 
     def sample(self, n):

tests/distributions/test_watson_negative_kappa_bingham.py

@@ -0,0 +1,48 @@
+import unittest
+
+import numpy.testing as npt
+import pyrecest.backend
+
+# pylint: disable=no-name-in-module,no-member
+from pyrecest.backend import array, linalg, ones
+from pyrecest.distributions import BinghamDistribution, WatsonDistribution
+
+
+@unittest.skipIf(
+    pyrecest.backend.__backend_name__ == "jax",
+    "Bingham conversion tests are not supported for this backend",
+)
+class TestWatsonNegativeKappaBingham(unittest.TestCase):
+    def test_negative_kappa_conversion_preserves_density(self):
+        mu = array([1.0, 0.0, 0.0, 0.0])
+        watson = WatsonDistribution(mu, -2.0)
+
+        bingham = watson.to_bingham()
+
+        self.assertIsInstance(bingham, BinghamDistribution)
+        npt.assert_allclose(bingham.Z, array([-2.0, 0.0, 0.0, 0.0]))
+        npt.assert_allclose(abs(float(bingham.M[:, 0] @ mu)), 1.0, atol=1e-7)
+
+        xs = array(
+            [
+                [1.0, 0.0, 0.0, 0.0],
+                [0.0, 1.0, 0.0, 0.0],
+                [1.0, 1.0, 0.0, 0.0],
+                [1.0, 1.0, 1.0, 1.0],
+            ]
+        )
+        xs = xs / linalg.norm(xs, axis=1).reshape((-1, 1))
+
+        npt.assert_allclose(watson.pdf(xs), bingham.pdf(xs), rtol=1e-6, atol=1e-8)
+
+    def test_negative_kappa_sampling_works_above_s2(self):
+        dist = WatsonDistribution(array([1.0, 0.0, 0.0, 0.0]), -2.0)
+
+        samples = dist.sample(4)
+
+        self.assertEqual(samples.shape, (4, 4))
+        npt.assert_allclose(linalg.norm(samples, axis=1), ones(4), atol=1e-6)
+
+
+if __name__ == "__main__":
+    unittest.main()