Merge pull request #18207 from dschmitz89/foldednormal_cdf

ENH: improve precision of folded normal distribution cdf

scipy/stats/_continuous_distns.py

@@ -2348,7 +2348,8 @@ def _pdf(self, x, c):
         return _norm_pdf(x + c) + _norm_pdf(x-c)
 
     def _cdf(self, x, c):
-        return _norm_cdf(x-c) + _norm_cdf(x+c) - 1.0
+        sqrt_two = np.sqrt(2)
+        return 0.5 * (sc.erf((x - c)/sqrt_two) + sc.erf((x + c)/sqrt_two))
 
     def _sf(self, x, c):
         return _norm_sf(x - c) + _norm_sf(x + c)

scipy/stats/tests/test_distributions.py

@@ -1122,6 +1122,19 @@ def test_entropy(self, c, ref):
         assert_allclose(stats.gompertz.entropy(c), ref, rtol=1e-14)
 
 
+class TestFoldNorm:
+
+    # reference values were computed with mpmath with 50 digits of precision
+    # from mpmath import mp
+    # mp.dps = 50
+    # mp.mpf(0.5) * (mp.erf((x - c)/mp.sqrt(2)) + mp.erf((x + c)/mp.sqrt(2)))
+
+    @pytest.mark.parametrize('x, c, ref', [(1e-4, 1e-8, 7.978845594730578e-05),
+                                           (1e-4, 1e-4, 7.97884555483635e-05)])
+    def test_cdf(self, x, c, ref):
+        assert_allclose(stats.foldnorm.cdf(x, c), ref, rtol=1e-15)
+
+
 class TestHalfNorm:
 
     # sfx is sf(x).  The values were computed with mpmath: