Fix double-counted hyperspherical moment Jacobian

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
 

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()