From 6923176471df0702f6b6821e14aa1f4abeb1bf44 Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Wed, 19 Aug 2026 17:47:18 +0800
Subject: [PATCH 1/2] Fix negative-kappa Watson Bingham conversion

---
 .../hypersphere_subset/watson_distribution.py | 22 +++++++++++++------
 1 file changed, 15 insertions(+), 7 deletions(-)

diff --git a/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py b/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py
index 9d2a650e2f..39b249c3a0 100644
--- a/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py
+++ b/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):

From 932e286ebae04772476a3f31ba4943a6f775ae21 Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Wed, 19 Aug 2026 17:47:40 +0800
Subject: [PATCH 2/2] Add negative-kappa Watson conversion regressions

---
 .../test_watson_negative_kappa_bingham.py     | 48 +++++++++++++++++++
 1 file changed, 48 insertions(+)
 create mode 100644 tests/distributions/test_watson_negative_kappa_bingham.py

diff --git a/tests/distributions/test_watson_negative_kappa_bingham.py b/tests/distributions/test_watson_negative_kappa_bingham.py
new file mode 100644
index 0000000000..ac6d0d2a23
--- /dev/null
+++ b/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()