From 63da5cc1504f066b31374027f637b4b021445d6b Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Fri, 25 Apr 2025 00:34:55 -0500 Subject: gh-132893: More accurate CDF computation (gh-132895) --- Lib/statistics.py | 9 +++++---- .../next/Library/2025-04-24-21-22-46.gh-issue-132893.KFuxZ2.rst | 1 + 2 files changed, 6 insertions(+), 4 deletions(-) create mode 100644 Misc/NEWS.d/next/Library/2025-04-24-21-22-46.gh-issue-132893.KFuxZ2.rst diff --git a/Lib/statistics.py b/Lib/statistics.py index e2b5926..9dcb20d 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 0000000..1ef6511 --- /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. -- cgit v0.12