From f318a088fef98538eb55062a01111e377f71b40c Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Sat, 29 Aug 2026 03:05:46 +0800
Subject: [PATCH 1/2] Add hyperspherical moment Jacobian regression

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

diff --git a/tests/distributions/test_hypersphere_moment_jacobian.py b/tests/distributions/test_hypersphere_moment_jacobian.py
new file mode 100644
index 000000000..351cf1586
--- /dev/null
+++ b/tests/distributions/test_hypersphere_moment_jacobian.py
@@ -0,0 +1,27 @@
+"""Regression tests for numerical hyperspherical moments."""
+
+import unittest
+
+import numpy as np
+import numpy.testing as npt
+import pyrecest.backend
+from pyrecest.distributions import HypersphericalUniformDistribution
+
+
+class HypersphereMomentJacobianTest(unittest.TestCase):
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ == "jax",
+        "Numerical hyperspherical integration is not supported on JAX.",
+    )
+    def test_uniform_s2_second_moment_uses_surface_jacobian_once(self):
+        """Uniform S2 has E[x x^T] = I/3 and therefore unit trace."""
+        dist = HypersphericalUniformDistribution(2)
+
+        moment = dist.moment_numerical()
+
+        npt.assert_allclose(moment, np.eye(3) / 3.0, atol=1e-6)
+        npt.assert_allclose(np.trace(moment), 1.0, atol=1e-6)
+
+
+if __name__ == "__main__":
+    unittest.main()

From 5507bb41689dcd5a32d8ddea22cb6988be07c3fa Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Sat, 29 Aug 2026 03:06:32 +0800
Subject: [PATCH 2/2] Apply hyperspherical surface Jacobian once in moments

---
 ...bstract_hypersphere_subset_distribution.py | 31 ++-----------------
 1 file changed, 3 insertions(+), 28 deletions(-)

diff --git a/src/pyrecest/distributions/hypersphere_subset/abstract_hypersphere_subset_distribution.py b/src/pyrecest/distributions/hypersphere_subset/abstract_hypersphere_subset_distribution.py
index 0ae1b395e..04eef9509 100644
--- a/src/pyrecest/distributions/hypersphere_subset/abstract_hypersphere_subset_distribution.py
+++ b/src/pyrecest/distributions/hypersphere_subset/abstract_hypersphere_subset_distribution.py
@@ -173,38 +173,13 @@ def f(points):
 
             return f
 
-        def g_gen(f_hypersph_coords, dim):
-            if dim == 1:
-
-                def g_1d(phi):
-                    return f_hypersph_coords(array(phi))
-
-                return g_1d
-            if dim == 2:
-
-                def g_2d(phi1, phi2):
-                    return f_hypersph_coords(array(phi1), array(phi2)) * sin(phi2)
-
-                return g_2d
-            if dim == 3:
-
-                def g_3d(phi1, phi2, phi3):
-                    return (
-                        f_hypersph_coords(array(phi1), array(phi2), array(phi3))
-                        * sin(phi2)
-                        * sin(phi3) ** 2
-                    )
-
-                return g_3d
-
-            raise ValueError("Dimension not supported.")
-
         for i in range(self.dim + 1):
             for j in range(self.dim + 1):
                 f_curr = f_gen(i, j)
                 fangles = self.__class__.gen_fun_hyperspherical_coords(f_curr, self.dim)
-                g_curr = g_gen(fangles, self.dim)
-                m[i, j] = self.__class__.integrate_fun_over_domain(g_curr, self.dim)
+                # integrate_fun_over_domain() already applies the hyperspherical
+                # surface-element Jacobian through integrate_fun_over_domain_part().
+                m[i, j] = self.__class__.integrate_fun_over_domain(fangles, self.dim)
 
         return m