Mailing List Archive

bpo-45876: Improve accuracy for stdev() and pstdev() in statistics (GH-29736)
https://github.com/python/cpython/commit/af9ee57b96cb872df6574e36027cc753417605f9
commit: af9ee57b96cb872df6574e36027cc753417605f9
branch: main
author: Raymond Hettinger <rhettinger@users.noreply.github.com>
committer: rhettinger <rhettinger@users.noreply.github.com>
date: 2021-11-26T22:54:50-07:00
summary:

bpo-45876: Improve accuracy for stdev() and pstdev() in statistics (GH-29736)

* Inlined code from variance functions

* Added helper functions for the float square root of a fraction

* Call helper functions

* Add blurb

* Fix over-specified test

* Add a test for the _sqrt_frac() helper function

* Increase the tested range

* Add type hints to the internal function.

* Fix test for correct rounding

* Simplify ??(n/m)? calculation

Co-authored-by: Mark Dickinson <dickinsm@gmail.com>

* Add comment and beef-up tests

* Test for zero denominator

* Add algorithmic references

* Add test for the _isqrt_frac_rto() helper function.

* Compute the 109 instead of hard-wiring it

* Stronger test for _isqrt_frac_rto()

* Bigger range

* Bigger range

* Replace float() call with int/int division to be parallel with the other code path.

* Factor out division. Update proof link. Remove internal type declaration

Co-authored-by: Mark Dickinson <dickinsm@gmail.com>

files:
A Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst
M Lib/statistics.py
M Lib/test/test_statistics.py

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 5c3f77df1549d..cf8eaa0a61e62 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -130,6 +130,7 @@
import math
import numbers
import random
+import sys

from fractions import Fraction
from decimal import Decimal
@@ -304,6 +305,27 @@ def _fail_neg(values, errmsg='negative value'):
raise StatisticsError(errmsg)
yield x

+def _isqrt_frac_rto(n: int, m: int) -> float:
+ """Square root of n/m, rounded to the nearest integer using round-to-odd."""
+ # Reference: https://www.lri.fr/~melquion/doc/05-imacs17_1-expose.pdf
+ a = math.isqrt(n // m)
+ return a | (a*a*m != n)
+
+# For 53 bit precision floats, the _sqrt_frac() shift is 109.
+_sqrt_shift: int = 2 * sys.float_info.mant_dig + 3
+
+def _sqrt_frac(n: int, m: int) -> float:
+ """Square root of n/m as a float, correctly rounded."""
+ # See principle and proof sketch at: https://bugs.python.org/msg407078
+ q = (n.bit_length() - m.bit_length() - _sqrt_shift) // 2
+ if q >= 0:
+ numerator = _isqrt_frac_rto(n, m << 2 * q) << q
+ denominator = 1
+ else:
+ numerator = _isqrt_frac_rto(n << -2 * q, m)
+ denominator = 1 << -q
+ return numerator / denominator # Convert to float
+

# === Measures of central tendency (averages) ===

@@ -837,14 +859,17 @@ def stdev(data, xbar=None):
1.0810874155219827

"""
- # Fixme: Despite the exact sum of squared deviations, some inaccuracy
- # remain because there are two rounding steps. The first occurs in
- # the _convert() step for variance(), the second occurs in math.sqrt().
- var = variance(data, xbar)
- try:
+ if iter(data) is data:
+ data = list(data)
+ n = len(data)
+ if n < 2:
+ raise StatisticsError('stdev requires at least two data points')
+ T, ss = _ss(data, xbar)
+ mss = ss / (n - 1)
+ if hasattr(T, 'sqrt'):
+ var = _convert(mss, T)
return var.sqrt()
- except AttributeError:
- return math.sqrt(var)
+ return _sqrt_frac(mss.numerator, mss.denominator)


def pstdev(data, mu=None):
@@ -856,14 +881,17 @@ def pstdev(data, mu=None):
0.986893273527251

"""
- # Fixme: Despite the exact sum of squared deviations, some inaccuracy
- # remain because there are two rounding steps. The first occurs in
- # the _convert() step for pvariance(), the second occurs in math.sqrt().
- var = pvariance(data, mu)
- try:
+ if iter(data) is data:
+ data = list(data)
+ n = len(data)
+ if n < 1:
+ raise StatisticsError('pstdev requires at least one data point')
+ T, ss = _ss(data, mu)
+ mss = ss / n
+ if hasattr(T, 'sqrt'):
+ var = _convert(mss, T)
return var.sqrt()
- except AttributeError:
- return math.sqrt(var)
+ return _sqrt_frac(mss.numerator, mss.denominator)


# === Statistics for relations between two inputs ===
diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index c0e427d9355f2..771a03e707ee0 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -9,13 +9,14 @@
import copy
import decimal
import doctest
+import itertools
import math
import pickle
import random
import sys
import unittest
from test import support
-from test.support import import_helper
+from test.support import import_helper, requires_IEEE_754

from decimal import Decimal
from fractions import Fraction
@@ -2161,6 +2162,66 @@ def test_center_not_at_mean(self):
self.assertEqual(self.func(data), 2.5)
self.assertEqual(self.func(data, mu=0.5), 6.5)

+class TestSqrtHelpers(unittest.TestCase):
+
+ def test_isqrt_frac_rto(self):
+ for n, m in itertools.product(range(100), range(1, 1000)):
+ r = statistics._isqrt_frac_rto(n, m)
+ self.assertIsInstance(r, int)
+ if r*r*m == n:
+ # Root is exact
+ continue
+ # Inexact, so the root should be odd
+ self.assertEqual(r&1, 1)
+ # Verify correct rounding
+ self.assertTrue(m * (r - 1)**2 < n < m * (r + 1)**2)
+
+ @requires_IEEE_754
+ def test_sqrt_frac(self):
+
+ def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
+ if not x:
+ return root == 0.0
+
+ # Extract adjacent representable floats
+ r_up: float = math.nextafter(root, math.inf)
+ r_down: float = math.nextafter(root, -math.inf)
+ assert r_down < root < r_up
+
+ # Convert to fractions for exact arithmetic
+ frac_root: Fraction = Fraction(root)
+ half_way_up: Fraction = (frac_root + Fraction(r_up)) / 2
+ half_way_down: Fraction = (frac_root + Fraction(r_down)) / 2
+
+ # Check a closed interval.
+ # Does not test for a midpoint rounding rule.
+ return half_way_down ** 2 <= x <= half_way_up ** 2
+
+ randrange = random.randrange
+
+ for i in range(60_000):
+ numerator: int = randrange(10 ** randrange(50))
+ denonimator: int = randrange(10 ** randrange(50)) + 1
+ with self.subTest(numerator=numerator, denonimator=denonimator):
+ x: Fraction = Fraction(numerator, denonimator)
+ root: float = statistics._sqrt_frac(numerator, denonimator)
+ self.assertTrue(is_root_correctly_rounded(x, root))
+
+ # Verify that corner cases and error handling match math.sqrt()
+ self.assertEqual(statistics._sqrt_frac(0, 1), 0.0)
+ with self.assertRaises(ValueError):
+ statistics._sqrt_frac(-1, 1)
+ with self.assertRaises(ValueError):
+ statistics._sqrt_frac(1, -1)
+
+ # Error handling for zero denominator matches that for Fraction(1, 0)
+ with self.assertRaises(ZeroDivisionError):
+ statistics._sqrt_frac(1, 0)
+
+ # The result is well defined if both inputs are negative
+ self.assertAlmostEqual(statistics._sqrt_frac(-2, -1), math.sqrt(2.0))
+
+
class TestStdev(VarianceStdevMixin, NumericTestCase):
# Tests for sample standard deviation.
def setUp(self):
@@ -2175,7 +2236,7 @@ def test_compare_to_variance(self):
# Test that stdev is, in fact, the square root of variance.
data = [random.uniform(-2, 9) for _ in range(1000)]
expected = math.sqrt(statistics.variance(data))
- self.assertEqual(self.func(data), expected)
+ self.assertAlmostEqual(self.func(data), expected)

def test_center_not_at_mean(self):
data = (1.0, 2.0)
diff --git a/Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst b/Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst
new file mode 100644
index 0000000000000..889ed6ce3ffb2
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst
@@ -0,0 +1,2 @@
+Improve the accuracy of stdev() and pstdev() in the statistics module. When
+the inputs are floats or fractions, the output is a correctly rounded float

_______________________________________________
Python-checkins mailing list
Python-checkins@python.org
https://mail.python.org/mailman/listinfo/python-checkins