8000 bpo-20499: Rounding error in statistics.pvariance (GH-28230) · python/cpython@4a5cccb · GitHub
[go: up one dir, main page]
Skip to content
< 10000 /div>

Commit 4a5cccb

Browse files
rhettingerrhettinger
authored
bpo-20499: Rounding error in statistics.pvariance (GH-28230)
1 parent f235dd0 commit 4a5cccb

File tree

3 files changed

+51
-56
lines changed

3 files changed

+51
-56
lines changed

‎Lib/statistics.py

Lines changed: 36 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -147,21 +147,17 @@ class StatisticsError(ValueError):
147147

148148
# === Private utilities ===
149149

150-
def _sum(data, start=0):
151-
"""_sum(data [, start]) -> (type, sum, count)
150+
def _sum(data):
151+
"""_sum(data) -> (type, sum, count)
152152
153153
Return a high-precision sum of the given numeric data as a fraction,
154154
together with the type to be converted to and the count of items.
155155
156-
If optional argument ``start`` is given, it is added to the total.
157-
If ``data`` is empty, ``start`` (defaulting to 0) is returned.
158-
159-
160156
Examples
161157
--------
162158
163-
>>> _sum([3, 2.25, 4.5, -0.5, 1.0], 0.75)
164-
(<class 'float'>, Fraction(11, 1), 5)
159+
>>> _sum([3, 2.25, 4.5, -0.5, 0.25])
160+
(<class 'float'>, Fraction(19, 2), 5)
165161
166162
Some sources of round-off error will be avoided:
167163
@@ -184,10 +180,9 @@ def _sum(data, start=0):
184180
allowed.
185181
"""
186182
count = 0
187-
n, d = _exact_ratio(start)
188-
partials = {d: n}
183+
partials = {}
189184
partials_get = partials.get
190-
T = _coerce(int, type(start))
185+
T = int
191186
for typ, values in groupby(data, type):
192187
T = _coerce(T, typ) # or raise TypeError
193188
for n, d in map(_exact_ratio, values):
@@ -200,8 +195,7 @@ def _sum(data, start=0):
200195
assert not _isfinite(total)
201196
else:
202197
# Sum all the partial sums using builtin sum.
203-
# FIXME is this faster if we sum them in order of the denominator?
204-
total = sum(Fraction(n, d) for d, n in sorted(partials.items()))
198+
total = sum(Fraction(n, d) for d, n in partials.items())
205199
return (T, total, count)
206200

207201

@@ -252,27 +246,19 @@ def _exact_ratio(x):
252246
x is expected to be an int, Fraction, Decimal or float.
253247
"""
254248
try:
255-
# Optimise the common case of floats. We expect that the most often
256-
# used numeric type will be builtin floats, so try to make this as
257-
# fast as possible.
258-
if type(x) is float or type(x) is Decimal:
259-
return x.as_integer_ratio()
260-
try:
261-
# x may be an int, Fraction, or Integral ABC.
262-
return (x.numerator, x.denominator)
263-
except AttributeError:
264-
try:
265-
# x may be a float or Decimal subclass.
266-
return x.as_integer_ratio()
267-
except AttributeError:
268-
# Just give up?
269-
pass
249+
return x.as_integer_ratio()
250+
except AttributeError:
251+
pass
270252
except (OverflowError, ValueError):
271253
# float NAN or INF.
272254
assert not _isfinite(x)
273255
return (x, None)
274-
msg = "can't convert type '{}' to numerator/denominator"
275-
raise TypeError(msg.format(type(x).__name__))
256+
try:
257+
# x may be an Integral ABC.
258+
return (x.numerator, x.denominator)
259+
except AttributeError:
260+
msg = f"can't convert type '{type(x).__name__}' to numerator/denominator"
261+
raise TypeError(msg)
276262

277263

278264
def _convert(value, T):
@@ -730,18 +716,20 @@ def _ss(data, c=None):
730716
if c is not None:
731717
T, total, count = _sum((d := x - c) * d for x in data)
732718
return (T, total)
733-
# Compute the mean accurate to within 1/2 ulp
734-
c = mean(data)
735-
# Initial computation for the sum of square deviations
736-
T, total, count = _sum((d := x - c) * d for x in data)
737-
# Correct any remaining inaccuracy in the mean c.
738-
# The following sum should mathematically equal zero,
739-
# but due to the final rounding of the mean, it may not.
740-
U, error, count2 = _sum((x - c) for x in data)
741-
assert count == count2
742-
correction = error * error / len(data)
743-
total -= correction
744-
assert not total < 0, 'negative sum of square deviations: %f' % total
719+
T, total, count = _sum(data)
720+
mean_n, mean_d = (total / count).as_integer_ratio()
721+
partials = Counter()
722+
for n, d in map(_exact_ratio, data):
723+
diff_n = n * mean_d - d * mean_n
724+
diff_d = d * mean_d
725+
partials[diff_d * diff_d] += diff_n * diff_n
726+
if None in partials:
727+
# The sum will be a NAN or INF. We can ignore all the finite
728+
# partials, and just look at this special one.
729+
total = partials[None]
730+
assert not _isfinite(total)
731+
else:
732+
total = sum(Fraction(n, d) for d, n in partials.items())
745733
return (T, total)
746734

747735

@@ -845,6 +833,9 @@ def stdev(data, xbar=None):
845833
1.0810874155219827
846834
847835
"""
836+
# Fixme: Despite the exact sum of squared deviations, some inaccuracy
837+
# remain because there are two rounding steps. The first occurs in
838+
# the _convert() step for variance(), the second occurs in math.sqrt().
848839
var = variance(data, xbar)
849840
try:
850841
return var.sqrt()
@@ -861,6 +852,9 @@ def pstdev(data, mu=None):
861852
0.986893273527251
862853
863854
"""
855+
# Fixme: Despite the exact sum of squared deviations, some inaccuracy
856+
# remain because there are two rounding steps. The first occurs in
857+
# the _convert() step for pvariance(), the second occurs in math.sqrt().
864858
var = pvariance(data, mu)
865859
try:
866860
return var.sqrt()

‎Lib/test/test_statistics.py

Lines changed: 14 additions & 14 deletions

‎Misc/NEWS.d/next/Library/2021-09-08-01-19-31.bpo-20499.tSxx8Y.rst

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1250,20 +1250,14 @@ def test_empty_data(self):
12501250
# Override test for empty data.
12511251
for data in ([], (), iter([])):
12521252
self.assertEqual(self.func(data), (int, Fraction(0), 0))
1253-
self.assertEqual(self.func(data, 23), (int, Fraction(23), 0))
1254-
self.assertEqual(self.func(data, 2.3), (float, Fraction(2.3), 0))
12551253

12561254
def test_ints(self):
12571255
self.assertEqual(self.func([1, 5, 3, -4, -8, 20, 42, 1]),
12581256
(int, Fraction(60), 8))
1259-
self.assertEqual(self.func([4, 2, 3, -8, 7], 1000),
1260-
(int, Fraction(1008), 5))
12611257

12621258
def test_floats(self):
12631259
self.assertEqual(self.func([0.25]*20),
12641260
(float, Fraction(5.0), 20))
1265-
self.assertEqual(self.func([0.125, 0.25, 0.5, 0.75], 1.5),
1266-
(float, Fraction(3.125), 4))
12671261

12681262
def test_fractions(self):
12691263
self.assertEqual(self.func([Fraction(1, 1000)]*500),
@@ -1284,14 +1278,6 @@ def test_compare_with_math_fsum(self):
12841278
data = [random.uniform(-100, 1000) for _ in range(1000)]
12851279
self.assertApproxEqual(float(self.func(data)[1]), math.fsum(data), rel=2e-16)
12861280

1287-
def test_start_argument(self):
1288-
# Test that the optional start argument works correctly.
1289-
data = [random.uniform(1, 1000) for _ in range(100)]
1290-
t = self.func(data)[1]
1291-
self.assertEqual(t+42, self.func(data, 42)[1])
1292-
self.assertEqual(t-23, self.func(data, -23)[1])
1293-
self.assertEqual(t+Fraction(1e20), self.func(data, 1e20)[1])
1294-
12951281
def test_strings_fail(self):
12961282
# Sum of strings should fail.
12971283
self.assertRaises(TypeError, self.func, [1, 2, 3], '999')
@@ -2101,6 +2087,13 @@ def test_decimals(self):
21012087
self.assertEqual(result, exact)
21022088
self.assertIsInstance(result, Decimal)
21032089

2090+
def test_accuracy_bug_20499(self):
2091+
data = [0, 0, 1]
2092+
exact = 2 / 9
2093+
result = self.func(data)
2094+
self.assertEqual(result, exact)
2095+
self.assertIsInstance(result, float)
2096+
21042097

21052098
class TestVariance(VarianceStdevMixin, NumericTestCase, UnivariateTypeMixin):
21062099
# Tests for sample variance.
@@ -2141,6 +2134,13 @@ def test_center_not_at_mean(self):
21412134
self.assertEqual(self.func(data), 0.5)
21422135
self.assertEqual(self.func(data, xbar=2.0), 1.0)
21432136

2137+
def test_accuracy_bug_20499(self):
2138+
data = [0, 0, 2]
2139+
exact = 4 / 3
2140+
result = self.func(data)
2141+
self.assertEqual(result, exact)
2142+
self.assertIsInstance(result, float)
2143+
21442144
class TestPStdev(VarianceStdevMixin, NumericTestCase):
21452145
# Tests for population standard deviation.
21462146
def setUp(self):
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
Improve the speed and accuracy of statistics.pvariance().

0 commit comments

Comments
 (0)
0