fix(distributions): improve lower tail numerical stability in Moyal logcdf (#8437)

pymc/distributions/continuous.py

@@ -4055,7 +4055,12 @@ def logp(value, mu, sigma):
 
     def logcdf(value, mu, sigma):
         scaled = (value - mu) / sigma
-        res = pt.log(pt.erfc(pt.exp(-scaled / 2) * (2**-0.5)))
+        u = pt.exp(-scaled / 2) * (2**-0.5)
+        res = pt.switch(
+            pt.gt(u, 1.0),
+            pt.log(pt.erfcx(u)) - pt.sqr(u),
+            pt.log(pt.erfc(u)),
+        )
         return check_parameters(
             res,
             sigma > 0,

pymc/gp/hsgp_approx.py

@@ -41,6 +41,11 @@ def set_boundary(X: TensorLike, c: numbers.Real | TensorLike) -> np.ndarray:
     S = (pt.max(X, axis=0) - pt.min(X, axis=0)) / 2.0
 
     L = (c * S).eval()  # eval() makes sure L is not changed with out-of-sample preds
+    if np.any(L <= 0):
+        raise ValueError(
+            "Boundary size L must be strictly positive for all dimensions. "
+            "Check that input data X has non-zero range (max - min > 0) for each feature column."
+        )
     return L
 
 

tests/distributions/test_continuous.py

@@ -956,6 +956,16 @@ def test_moyal_logcdf(self):
         if pytensor.config.floatX == "float32":
             raise Exception("Flaky test: It passed this time, but XPASS is not allowed.")
 
+    def test_moyal_logcdf_lower_tail(self):
+        rv = pm.Moyal.dist(mu=0.0, sigma=1.0)
+        xs = np.array([-10.0, -20.0, -50.0])
+        res = pm.logcdf(rv, xs).eval()
+        # F(x) = 2 * Phi(-exp(-x / 2)); avoid taking log of an underflowed CDF.
+        expected = np.log(2.0) + sp.log_ndtr(-np.exp(-xs / 2.0))
+        assert np.isfinite(expected).all()
+        assert np.isfinite(res).all()
+        np.testing.assert_allclose(res, expected, rtol=1e-5)
+
     def test_moyal_icdf(self):
         check_icdf(
             pm.Moyal,

tests/gp/test_hsgp_approx.py

@@ -121,6 +121,11 @@ def test_set_boundaries_3d(self, X2):
         L = pm.gp.hsgp_approx.set_boundary(X2s, c=2)
         assert np.all(L == 10)
 
+    def test_set_boundary_zero_range_raises(self):
+        X = np.ones((10, 2))
+        with pytest.raises(ValueError, match="Boundary size L must be strictly positive"):
+            pm.gp.hsgp_approx.set_boundary(X, c=1.5)
+
     def test_mean_invariance(self):
         X = np.linspace(0, 10, 100)[:, None]
         original_center = (np.max(X, axis=0) - np.min(X, axis=0)) / 2