Repository: dask/dask

Issue:
`quantile` calculation using dask default method gives incorrect results
This is similar to #731 but reveals a broader issue and thus I would like to create a new issue for it. I'm using `dask==2.24.0` and `pandas==1.0.5`
```
import pandas as pd
import dask.dataframe as dd
x = pd.DataFrame({'a': [1,2,3]})
xd = dd.from_pandas(x, npartitions=3)
print(x['a'].quantile(0.5))
print(xd['a'].quantile(0.5).compute())
```

This difference persists as I try all natural number column `'a'` with length greater than 2 (i.e. `'a'` being natural numbers 1,...,N where N>2).

```
for i in range(1, 100):
    x = pd.DataFrame({'a': list(range(i))})
    xd = dd.from_pandas(x, npartitions=i+1)
    pd_quantile = x['a'].quantile(0.5)
    dd_quantile = xd['a'].quantile(0.5).compute()
    if pd_quantile == dd_quantile:
        print('same')
    else:
        print('different')
```

Using `method='tdigest'` does solve the issue, but the default algorithm shouldn't be giving a wrong answer either..

Discussion:
mrocklin: cc @jcrist if you have time
jcrist: From the issue https://github.com/dask/dask/issues/6551#issuecomment-681964573 it looks like this doesn't fully resolve the issue. I'm skeptical of work on this front - the dask method is custom and lacks the statistical rigor users might expect. I'm skeptical of one-off fixes that improve things for certain use cases, as without statistical rigor the results are likely made worse for other users. Perhaps we should switch to using `tdigest` as the default, since that's much more robust?
eriknw: I agree on most counts.  We should use `tdigest` if possible.  I am skeptical of quick fixes.  The current method was never meant for generic percentile calculations and lacks the statistical rigor to do so.  Additional work (beyond what is in this PR IMHO) to improve or investigate the current approach for quantiles is probably not worth the effort.

As the person who wrote the code for *partition quantiles* (i.e., determining how to partition dask dataframes, which it does an excellent job at), it always bothered me how it was adapted for generic percentiles, and, in particular, how it doesn't include the 100% percentile.  This nags me, because I have knowledge of the implementation.  It should include 100%, and including 0% should also help at some boundary cases.

So, as somebody who understands the current implementation, I trust that this PR generally makes things better.  I'm sorry for not making it years ago.

This PR also *substantially improves* #6551 to the point where I think it's resolved within expectation.  I think the additional edge cases that were found should be expected, because these are *approximate* quantiles, and the default implementation doesn't have enough information to correctly gather and merge the partition summaries.

@jcrist shall I update the PR to make `tdigest` default if available, or do you still disagree in principle of making any changes to the current implementation?  The current implementation will still be used for non-numeric data, so I think improvements to the current implementation should be accepted if they are available.
jcrist: I'm fine with improving the current implementation, but think we should start a discussion about no longer defaulting to it. I'll open a new issue for that. The fix you have provided here and reasoning about it seems good to me, so I'm happy to merge (looks like the new black release has broke things here too though, black changed how it handles trailing commas).
jcrist: See #6566

Patch:
diff --git a/dask/array/percentile.py b/dask/array/percentile.py
index ef4e9c09f..e35c6f7a0 100644
--- a/dask/array/percentile.py
+++ b/dask/array/percentile.py
@@ -126,10 +126,12 @@ def percentile(a, q, interpolation="linear", method="default"):
 
     # Otherwise use the custom percentile algorithm
     else:
-
+        # Add 0 and 100 during calculation for more robust behavior (hopefully)
+        calc_q = np.pad(q, 1, mode="constant")
+        calc_q[-1] = 100
         name = "percentile_chunk-" + token
         dsk = dict(
-            ((name, i), (_percentile, key, q, interpolation))
+            ((name, i), (_percentile, key, calc_q, interpolation))
             for i, key in enumerate(a.__dask_keys__())
         )
 
@@ -138,7 +140,7 @@ def percentile(a, q, interpolation="linear", method="default"):
             (name2, 0): (
                 merge_percentiles,
                 q,
-                [q] * len(a.chunks[0]),
+                [calc_q] * len(a.chunks[0]),
                 sorted(dsk),
                 interpolation,
             )
diff --git a/dask/dataframe/core.py b/dask/dataframe/core.py
index b8a134c5f..18667b6ce 100644
--- a/dask/dataframe/core.py
+++ b/dask/dataframe/core.py
@@ -5457,16 +5457,19 @@ def quantile(df, q, method="default"):
 
         from dask.array.percentile import _percentile, merge_percentiles
 
+        # Add 0 and 100 during calculation for more robust behavior (hopefully)
+        calc_qs = np.pad(qs, 1, mode="constant")
+        calc_qs[-1] = 100
         name = "quantiles-1-" + token
         val_dsk = {
-            (name, i): (_percentile, (getattr, key, "values"), qs)
+            (name, i): (_percentile, (getattr, key, "values"), calc_qs)
             for i, key in enumerate(df.__dask_keys__())
         }
 
         name2 = "quantiles-2-" + token
         merge_dsk = {
             (name2, 0): finalize_tsk(
-                (merge_percentiles, qs, [qs] * df.npartitions, sorted(val_dsk))
+                (merge_percentiles, qs, [calc_qs] * df.npartitions, sorted(val_dsk))
             )
         }
     dsk = merge(val_dsk, merge_dsk)


Test patch:
diff --git a/dask/array/tests/test_array_function.py b/dask/array/tests/test_array_function.py
index a76530bce..f016a6a53 100644
--- a/dask/array/tests/test_array_function.py
+++ b/dask/array/tests/test_array_function.py
@@ -179,7 +179,13 @@ def test_non_existent_func():
 
 @pytest.mark.skipif(missing_arrfunc_cond, reason=missing_arrfunc_reason)
 @pytest.mark.parametrize(
-    "func", [np.equal, np.matmul, np.dot, lambda x, y: np.stack([x, y])]
+    "func",
+    [
+        np.equal,
+        np.matmul,
+        np.dot,
+        lambda x, y: np.stack([x, y]),
+    ],
 )
 @pytest.mark.parametrize(
     "arr_upcast, arr_downcast",
diff --git a/dask/array/tests/test_percentiles.py b/dask/array/tests/test_percentiles.py
index c8669d6c2..fc118b426 100644
--- a/dask/array/tests/test_percentiles.py
+++ b/dask/array/tests/test_percentiles.py
@@ -80,6 +80,15 @@ def test_percentiles_with_empty_arrays(method):
     )
 
 
+@percentile_methods
+def test_percentiles_with_empty_q(method):
+    x = da.ones(10, chunks=((5, 0, 5),))
+    assert_eq(
+        da.percentile(x, [], method=method),
+        np.array([], dtype=x.dtype),
+    )
+
+
 @percentile_methods
 @pytest.mark.parametrize("q", [5, 5.0, np.int64(5), np.float64(5)])
 def test_percentiles_with_scaler_percentile(method, q):
diff --git a/dask/dataframe/tests/test_dataframe.py b/dask/dataframe/tests/test_dataframe.py
index 8b2aca252..9287c4455 100644
--- a/dask/dataframe/tests/test_dataframe.py
+++ b/dask/dataframe/tests/test_dataframe.py
@@ -1190,7 +1190,7 @@ def test_nbytes():
 
 @pytest.mark.parametrize(
     "method,expected",
-    [("tdigest", (0.35, 3.80, 2.5, 6.5, 2.0)), ("dask", (0.0, 5.4, 1.2, 7.8, 5.0))],
+    [("tdigest", (0.35, 3.80, 2.5, 6.5, 2.0)), ("dask", (0.0, 4.0, 1.2, 6.2, 2.0))],
 )
 def test_quantile(method, expected):
     if method == "tdigest":
@@ -1278,9 +1278,9 @@ def test_empty_quantile(method):
         (
             "dask",
             (
-                pd.Series([16.5, 36.5, 26.5], index=["A", "X", "B"]),
+                pd.Series([7.0, 27.0, 17.0], index=["A", "X", "B"]),
                 pd.DataFrame(
-                    [[1.50, 21.50, 11.50], [17.75, 37.75, 27.75]],
+                    [[1.50, 21.50, 11.50], [14.0, 34.0, 24.0]],
                     index=[0.25, 0.75],
                     columns=["A", "X", "B"],
                 ),
@@ -1355,6 +1355,14 @@ def test_quantile_for_possibly_unsorted_q():
     assert_eq(r, 25.0)
 
 
+def test_quantile_tiny_partitions():
+    """ See https://github.com/dask/dask/issues/6551 """
+    df = pd.DataFrame({"a": [1, 2, 3]})
+    ddf = dd.from_pandas(df, npartitions=3)
+    r = ddf["a"].quantile(0.5).compute()
+    assert r == 2
+
+
 def test_index():
     assert_eq(d.index, full.index)