Skip to content

statistics.quantiles(exclusive) returns unsorted cut points for duplicate float-inexact data #150318

Open
Open
statistics.quantiles(exclusive) returns unsorted cut points for duplicate float-inexact data#150318
Assignees
rhettinger
Labels
stdlibStandard Library Python modules in the Lib/ directorytype-bugAn unexpected behavior, bug, or error

Description

@zhangjiashuo-cs
zhangjiashuo-cs
opened on May 23, 2026

Bug description

statistics.quantiles(data, n=N, method='exclusive') with N ≥ 3 returns cut points that are not monotonically non-decreasing when data contains two equal but float-inexact values. The first and last cut points compute one ULP below the data value, while interior cut points compute exactly the data value, so the returned list is unsorted.

import statistics

# Two copies of float(pi)
data = [3.141592653589793, 3.141592653589793]
result = statistics.quantiles(data, n=9, method='exclusive')

print(result)
# [3.1415926535897927, 3.141592653589793, 3.141592653589793,
#  3.141592653589793, 3.141592653589793, 3.141592653589793,
#  3.141592653589793, 3.1415926535897927]

print(result == sorted(result))
# False  <-- cut points 0 and 7 are < cut points 1..6

Reproducer is one-liner:

>>> import statistics
>>> r = statistics.quantiles([3.141592653589793]*2, n=9, method='exclusive')
>>> r[0] <= r[1]
False

Why this is a bug

The quantiles documentation says it "Divide(s) data into n continuous intervals with equal probability. Returns a list of n - 1 cut points separating the intervals." Cut points that separate intervals must be monotonically non-decreasing — otherwise the "intervals" they delimit overlap or invert. Several downstream uses (e.g. bisect.bisect_left(cut_points, x) to bucket a value) rely silently on this invariant and will return wrong bucket indices when it's violated.

The 'exclusive' method computes each cut point with the formula:

# Lib/statistics.py
m = ld - 1
result = []
for i in range(1, n):
    j = i * m // n
    delta = i * m - j * n
    interpolated = (data[j-1] * (n - delta) + data[j] * delta) / n
    result.append(interpolated)

When data[j-1] == data[j] == x exactly, interpolated = (x*(n-delta) + x*delta)/n should be x mathematically, but for non-power-of-two delta/n ratios the float operations x*(n-delta) + x*delta do not always round to n*x, yielding a result one ULP off. The error happens asymmetrically across i, so adjacent cut points differ by 1 ULP and the list ends up unsorted.

'inclusive' is not affected because its formula is structured to be exact when data are identical.

Suggested fix

Sort the result before returning (cheap, O(n log n) on n − 1 elements which is already trivial):

--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -X,Y +X,Y @@ def quantiles(data, *, n=4, method='exclusive'):
         for i in range(1, n):
             j = i * m // n             # rescale i to m/n
             delta = i*m - j*n          # exact remainder
             interpolated = (data[j-1] * (n - delta) + data[j] * delta) / n
             result.append(interpolated)
-        return result
+        # The interpolation above can produce float-rounding artifacts that
+        # leave neighbouring identical-data cut points one ULP out of order
+        # (e.g. data=[pi, pi], n=9).  Restore the documented monotonic-
+        # non-decreasing invariant before returning.
+        result.sort()
+        return result

This preserves the existing values exactly when they are already monotonic and only reorders the few corner cases where float rounding broke the invariant. A regression test:

import statistics
def test_quantiles_monotonic_with_duplicates():
    for x in [3.141592653589793, 1/3, 0.1, 1e300, 1e-300]:
        for n in range(2, 20):
            for method in ('exclusive', 'inclusive'):
                r = statistics.quantiles([x]*2, n=n, method=method)
                assert r == sorted(r), (x, n, method, r)

CPython versions tested on

3.9, 3.12, 3.14

Operating systems tested on

macOS

Linked PRs

Activity

  1. aisk
    added
    type-bugAn unexpected behavior, bug, or error
    stdlibStandard Library Python modules in the Lib/ directory
    on May 23, 2026

  2. edvilme commented on May 23, 2026

    @edvilme
    edvilme
    on May 23, 2026
    Contributor

    Hi. If you're not planning on opening a PR, I would love to work on this 😁

  3. JelleZijlstra commented on May 23, 2026

    @JelleZijlstra
    JelleZijlstra
    on May 23, 2026
    Member

    Is sorting really the right answer here?

  4. rhettinger
    self-assigned this
    on May 23, 2026

  5. rhettinger commented on May 23, 2026

    @rhettinger
    rhettinger
    on May 23, 2026
    Contributor

    do not always round to n*x, yielding a result one ULP off.

    In practice, I don't think this would ever matter.

  6. edvilme commented on May 23, 2026

    @edvilme
    edvilme
    on May 23, 2026
    Contributor

    Is sorting really the right answer here?

    I think it would be better to catch the condition here, otherwise it would sort unnecesarily

  7. skirpichev commented on May 24, 2026

    @skirpichev
    skirpichev
    on May 24, 2026
    Member

    Is sorting really the right answer here?

    I think it's a bad answer.

    I think it would be better to catch the condition here

    How about using different (equivalent for reals) expression:

    interpolated = data[j-1] + (data[j] - data[j-1])*delta/n

    ?

  8. lbaumbauer commented on Jul 30, 2026

    @lbaumbauer
    lbaumbauer
    on Jul 30, 2026

    'inclusive' is not affected

    I am getting the same error with the method inclusive.

        assert r == sorted(r), (x, n, method, r)
           ^^^^^^^^^^^^^^
AssertionError: (3.141592653589793, 13, 'inclusive', [3.1415926535897936, 3.1415926535897927, 3.1415926535897936, 3.1415926535897936, 3.1415926535897936, 3.1415926535897936, 3.1415926535897936, 3.1415926535897936, 3.1415926535897936, 3.1415926535897936, 3.1415926535897927, 3.1415926535897936])

    Inclusive Code:

        m = ld - 1
        result = []
        for i in range(1, n):
            j, delta = divmod(i * m, n)
            interpolated = (data[j] * (n - delta) + data[j + 1] * delta) / n
            result.append(interpolated)
        return result

    Test case:

     def test_quantiles_monotonic_with_duplicates():
    for x in [3.141592653589793, 1/3, 0.1, 1e300, 1e-300]:
        for n in range(2, 20):
            for method in ('exclusive', 'inclusive'):
                r = quantiles([x]*2, n=n, method=method)
                assert r == sorted(r), (x, n, method, r)

test_quantiles_monotonic_with_duplicates()

    Fixed code (not 100% the same code as the suggested fix in the mentioned commit, just a demo with a single if clause added, does not handle delta issues or float)

        if method == 'inclusive':
        m = ld - 1
        result = []
        for i in range(1, n):
            j, delta = divmod(i * m, n)
            if data[j] == data[j + 1]:
                interpolated = data[j]
            else:
                interpolated = (data[j] * (n - delta) + data[j + 1] * delta) / n
            result.append(interpolated)
        return result

    Edit: fixed j, j+1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

Labels

stdlibStandard Library Python modules in the Lib/ directorytype-bugAn unexpected behavior, bug, or error

Projects

No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions