Merge pull request numpy#20331 from bzah/fix/quantile-default-lerp-method

MAINT: Update quantile default lerp method

Author: seberg

numpy/lib/function_base.py

@@ -67,7 +67,7 @@
 # fix_gamma : Callable
 #   A function used for discret methods to force the index to a specific value.
 _QuantileInterpolation = dict(
-    # --- HYNDMAN and FAN methods
+    # --- HYNDMAN AND FAN METHODS
     # Discrete methods
     inverted_cdf=dict(
         get_virtual_index=lambda n, quantiles: _inverted_cdf(n, quantiles),
@@ -102,10 +102,12 @@
         _compute_virtual_index(n, quantiles, 0, 0),
         fix_gamma=lambda gamma, _: gamma,
     ),
-    # Default value
+    # Default method.
+    # To avoid some rounding issues, `(n-1) * quantiles` is preferred to
+    # `_compute_virtual_index(n, quantiles, 1, 1)`.
+    # They are mathematically equivalent.
     linear=dict(
-        get_virtual_index=lambda n, quantiles:
-        _compute_virtual_index(n, quantiles, 1, 1),
+        get_virtual_index=lambda n, quantiles: (n - 1) * quantiles,
         fix_gamma=lambda gamma, _: gamma,
     ),
     median_unbiased=dict(
@@ -118,7 +120,7 @@
         _compute_virtual_index(n, quantiles, 3 / 8.0, 3 / 8.0),
         fix_gamma=lambda gamma, _: gamma,
     ),
-    # --- OTHER METHODS fixme add deprecated ?
+    # --- OTHER METHODS
     lower=dict(
         get_virtual_index=lambda n, quantiles: np.floor(
             (n - 1) * quantiles).astype(np.intp),

numpy/lib/tests/test_function_base.py

@@ -3356,6 +3356,14 @@ def test_nan_q(self):
 class TestQuantile:
     # most of this is already tested by TestPercentile
 
+    def test_max_ulp(self):
+        x = [0.0, 0.2, 0.4]
+        a = np.quantile(x, 0.45)
+        # The default linear method would result in 0 + 0.2 * (0.45/2) = 0.18.
+        # 0.18 is not exactly representable and the formula leads to a 1 ULP
+        # different result. Ensure it is this exact within 1 ULP, see gh-20331.
+        np.testing.assert_array_max_ulp(a, 0.18, maxulp=1)
+
     def test_basic(self):
         x = np.arange(8) * 0.5
         assert_equal(np.quantile(x, 0), 0.)