From: Raymond Hettinger <rhettinger@users.noreply.github.com>
Date: Fri, 25 Apr 2025 05:34:55 +0000 (-0500)
Subject: gh-132893: More accurate CDF computation (gh-132895)
X-Git-Tag: v3.14.0b1~306
X-Git-Url: http://git.ipfire.org/gitweb.cgi?a=commitdiff_plain;h=63da5cc1504f066b31374027f637b4b021445d6b;p=thirdparty%2FPython%2Fcpython.git

gh-132893: More accurate CDF computation (gh-132895)
---

diff --git a/Lib/statistics.py b/Lib/statistics.py
index e2b59267f04f..9dcb20d2a485 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -138,7 +138,7 @@ from fractions import Fraction
 from decimal import Decimal
 from itertools import count, groupby, repeat
 from bisect import bisect_left, bisect_right
-from math import hypot, sqrt, fabs, exp, erf, tau, log, fsum, sumprod
+from math import hypot, sqrt, fabs, exp, erfc, tau, log, fsum, sumprod
 from math import isfinite, isinf, pi, cos, sin, tan, cosh, asin, atan, acos
 from functools import reduce
 from operator import itemgetter
@@ -811,8 +811,9 @@ def register(*kernels):
 def normal_kernel():
     sqrt2pi = sqrt(2 * pi)
     sqrt2 = sqrt(2)
+    neg_sqrt2 = -sqrt2
     pdf = lambda t: exp(-1/2 * t * t) / sqrt2pi
-    cdf = lambda t: 1/2 * (1.0 + erf(t / sqrt2))
+    cdf = lambda t: 1/2 * erfc(t / neg_sqrt2)
     invcdf = lambda t: _normal_dist_inv_cdf(t, 0.0, 1.0)
     support = None
     return pdf, cdf, invcdf, support
@@ -1257,7 +1258,7 @@ class NormalDist:
         "Cumulative distribution function.  P(X <= x)"
         if not self._sigma:
             raise StatisticsError('cdf() not defined when sigma is zero')
-        return 0.5 * (1.0 + erf((x - self._mu) / (self._sigma * _SQRT2)))
+        return 0.5 * erfc((self._mu - x) / (self._sigma * _SQRT2))
 
     def inv_cdf(self, p):
         """Inverse cumulative distribution function.  x : P(X <= x) = p
@@ -1311,7 +1312,7 @@ class NormalDist:
         dv = Y_var - X_var
         dm = fabs(Y._mu - X._mu)
         if not dv:
-            return 1.0 - erf(dm / (2.0 * X._sigma * _SQRT2))
+            return erfc(dm / (2.0 * X._sigma * _SQRT2))
         a = X._mu * Y_var - Y._mu * X_var
         b = X._sigma * Y._sigma * sqrt(dm * dm + dv * log(Y_var / X_var))
         x1 = (a + b) / dv
diff --git a/Misc/NEWS.d/next/Library/2025-04-24-21-22-46.gh-issue-132893.KFuxZ2.rst b/Misc/NEWS.d/next/Library/2025-04-24-21-22-46.gh-issue-132893.KFuxZ2.rst
new file mode 100644
index 000000000000..1ef651188e04
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2025-04-24-21-22-46.gh-issue-132893.KFuxZ2.rst
@@ -0,0 +1 @@
+Improved the accuracy of ``statistics.NormalDist.cdf`` for negative inputs.