Improve accuracy of stdev functions in statistics

[rhettinger]

BPO	45876
Nosy	@tim-one, @rhettinger, @mdickinson, @stevendaprano
PRs	bpo-45876: Improve accuracy for stdev() and pstdev() in statistics #29736 bpo-45876: Correctly rounded stdev() and pstdev() for the Decimal case #29828 bpo-45876: Have stdev() also use decimal specific square root. #29869

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = 'https://github.com/rhettinger'
closed_at = <Date 2021-11-27.05:55:16.330>
created_at = <Date 2021-11-23.08:21:50.548>
labels = ['type-feature', 'library', '3.11']
title = 'Improve accuracy of stdev functions in statistics'
updated_at = <Date 2021-12-01.01:26:07.155>
user = 'https://github.com/rhettinger'

bugs.python.org fields:

activity = <Date 2021-12-01.01:26:07.155>
actor = 'rhettinger'
assignee = 'rhettinger'
closed = True
closed_date = <Date 2021-11-27.05:55:16.330>
closer = 'rhettinger'
components = ['Library (Lib)']
creation = <Date 2021-11-23.08:21:50.548>
creator = 'rhettinger'
dependencies = []
files = []
hgrepos = []
issue_num = 45876
keywords = ['patch']
message_count = 36.0
messages = ['406824', '406830', '406843', '406850', '406859', '406866', '406870', '406875', '406885', '406886', '406889', '406896', '406897', '406911', '406940', '406971', '407024', '407025', '407031', '407037', '407060', '407063', '407078', '407083', '407095', '407096', '407097', '407098', '407100', '407125', '407129', '407133', '407134', '407135', '407415', '407420']
nosy_count = 4.0
nosy_names = ['tim.peters', 'rhettinger', 'mark.dickinson', 'steven.daprano']
pr_nums = ['29736', '29828', '29869']
priority = 'normal'
resolution = 'fixed'
stage = 'resolved'
status = 'closed'
superseder = None
type = 'enhancement'
url = 'https://bugs.python.org/issue45876'
versions = ['Python 3.11']

[rhettinger]

The standard deviation computation in the statistics module is still subject to error even though the mean and sum of square differences are computed exactly using fractions.

The problem is that the exact fraction gets rounded to a float before going into math.sqrt() which also rounds.

It would be better to compute a high precision square root directly from the exact fraction rather than from a fraction rounded to a float.

Here are two ways to do it. With Cpython, the second way is twice as fast. With PyPy3, the speed is about the same.

def frac_sqrt(x: Fraction) -> float:
    'Return sqrt to greater precision than math.sqrt()'
    # Needed because a correctly rounded square root of
    # a correctly rounded input can still be off by 1 ulp.
    # Here we avoid the initial rounding and work directly
    # will the exact fractional input.  The square root
    # is first approximated with math.sqrt() and then
    # refined with a divide-and-average step. Since the
    # Newton-Raphson algorithm has quadratic convergence,
    # one refinement step gives excellent accuracy.
    a = Fraction(math.sqrt(x))
    return float((x / a + a) / 2)


def deci_sqrt(x: Fraction) -> float:
    ratio = Decimal(x.numerator) / Decimal(x.denominator)
    return float(ratio.sqrt())

[mdickinson]

I'm not sure this is worth worrying about. We already have a very tight error bound on the result: if x is a (positive) fraction and y is the closest float to x, (and assuming IEEE 754 binary64, round-ties-to-even, no overflow or underflow, etc.) then math.sqrt(y) will be in error by strictly less than 1 ulp from the true value √x, so we're already faithfully rounded. (And in particular, if the std. dev. is exactly representable as a float, this guarantees that we'll get that standard deviation exactly.)

[mdickinson]

One thought: would the deci_sqrt approach help with value ranges where the values are well within float limits, but the squares of the values are not? E.g., on my machine, I currently get errors for both of the following:

>> xs = [random.normalvariate(0.0, 1e200) for _ in range(10**6)]
>> statistics.stdev(xs)

>> xs = [random.normalvariate(0.0, 1e-200) for _ in range(10**6)]
>> statistics.stdev(xs)

It's hard to imagine that there are too many use-cases for values of this size, but it still feels a bit odd to be constrained to only half of the dynamic range of float.

[rhettinger]

I'm not sure this is worth worrying about ...

Instead of writing simple, mostly accurate code with math.fsum(), these functions have already applied labor intensive measures to get an exact mean and exact sum of square differences expressed in fractions. Then at the final step before the square root, it prematurely rounds to a float. This is easy to fix and has a single step cost comparable to that already paid for each input datum.

In a sports analogy, we've run the ball almost the full length of the field and then failed to put the ball over the goal line.

Part of the module's value proposition is that it strives to be more accurate than obvious implementations. The mean, variance, and pvariance function are correctly rounded. In this regard, the stdev and pstdev functions are deficient, but they could easily be made to be almost always correctly rounded.

One thought: would the deci_sqrt approach help with value ranges
where the values are well within float limits, but the squares
of the values are not?

Yes, the Emin and Emax for the default context is already almost big enough:

    Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, 
    Emax=999999, capitals=1, clamp=0, flags=[], 
    traps=[InvalidOperation, DivisionByZero, Overflow])

We could bump that up to fully envelop operations on the sum of squares:

    Context(Emin=-9_999_999, Emax=9_999_999, prec=50)

[mdickinson]

we've run the ball almost the full length of the field and then failed to put the ball over the goal line

But if we only go from "faithfully rounded" to "almost always correctly rounded", it seems to me that we're still a couple of millimetres away from that goal line. It wouldn't be hard to go for _always_ correctly rounded and actually get it over.

Yes, the Emin and Emax for the default context is already almost big enough

I'm confused: big enough for what? I was thinking of the use-case where the inputs are all floats, in which case an Emax of 999 and an Emin of -999 would already be more than big enough.

[rhettinger]

It wouldn't be hard to go for _always_ correctly rounded
and actually get it over.

Yes, that would be the right thing to do.

Does the technique you had in mind involve testing 1 ulp up or down to see whether its square is closer to the input?

an Emax of 999 and an Emin of -999 would already be
more than big enough.

Right. In haste, I confused max_exp=1024 with max_10_exp=308.

Am still thinking that the precision needs to be bumped up a bit the 28 place default.

[mdickinson]

Does the technique you had in mind involve testing 1 ulp up or down to see whether its square is closer to the input?

Kinda sorta. Below is some code: it's essentially just pure integer operations, with a last minute conversion to float (implicit in the division in the case of the second branch). And it would need to be better tested, documented, and double-checked to be viable.

def isqrt_rto(n):
    """
    Square root of n, rounded to the nearest integer using round-to-odd.
    """
    a = math.isqrt(n)
    return a | (a*a != n)


def isqrt_frac_rto(n, m):
    """
    Square root of n/m, rounded to the nearest integer using round-to-odd.
    """
    quotient, remainder = divmod(isqrt_rto(4*n*m), 2*m)
    return quotient | bool(remainder)


def sqrt_frac(n, m):
    """
    Square root of n/m as a float, correctly rounded.
    """
    quantum = (n.bit_length() - m.bit_length() - 1) // 2 - 54
    if quantum >= 0:
        return float(isqrt_frac_rto(n, m << 2 * quantum) << quantum)
    else:
        return isqrt_frac_rto(n << -2 * quantum, m) / (1 << -quantum)

[mdickinson]

Here's the float-and-Fraction-based code that I'm using to compare the integer-based code against:

def sqrt_frac2(n, m):
    """
    Square root of n/m as a float, correctly rounded.
    """
    f = fractions.Fraction(n, m)

    # First approximation.
    x = math.sqrt(n / m)

    # Use the approximation to find a pair of floats bracketing the actual sqrt
    if fractions.Fraction(x)**2 >= f:
        x_lo, x_hi = math.nextafter(x, 0.0), x
    else:
        x_lo, x_hi = x, math.nextafter(x, math.inf)

# Check the bracketing. If math.sqrt is correctly rounded (as it will be on a
# typical machine), then the assert can't fail. But we can't rely on math.sqrt being
# correctly rounded in general, so would need some fallback.
fx_lo, fx_hi = fractions.Fraction(x_lo), fractions.Fraction(x_hi)
assert fx_lo**2 <= f <= fx_hi**2

    # Compare true square root with the value halfway between the two floats.
    mid = (fx_lo + fx_hi) / 2
    if mid**2 < f:
        return x_hi
    elif mid**2 > f:
        return x_lo
    else:
        # Tricky case: mid**2 == f, so we need to choose the "even" endpoint.
        # Cheap trick: the addition in 0.5 * (x_lo + x_hi) will round to even.
        return 0.5 * (x_lo + x_hi)

[rhettinger]

Should the last line of sqrt_frac() be wrapped with float()?

[mdickinson]

Should the last line of sqrt_frac() be wrapped with float()?

It's already a float - it's the result of an int / int division.

[mdickinson]

Hmm. isqrt_frac_rto is unnecessarily complicated. Here's a more streamlined version of the code.

import math

def isqrt_frac_rto(n, m):
    """
    Square root of n/m, rounded to the nearest integer using round-to-odd.
    """
    a = math.isqrt(n*m) // m
    return a | (a*a*m != n)

def sqrt_frac(n, m):
    """
    Square root of n/m as a float, correctly rounded.
    """
    q = (n.bit_length() - m.bit_length() - 109) // 2
    if q >= 0:
        return float(isqrt_frac_rto(n, m << 2 * q) << q)
    else:
        return isqrt_frac_rto(n << -2 * q, m) / (1 << -q)