Merge pull request numpy#26769 from aherbert/fix-closest-observation

BUG: Quantile closest_observation to round to nearest even order

Author: charris

doc/release/upcoming_changes/26656.improvement.rst

@@ -0,0 +1,5 @@
+`np.quantile` with method ``closest_observation`` chooses nearest even order statistic
+--------------------------------------------------------------------------------------
+This changes the definition of nearest for border cases from the nearest odd
+order statistic to nearest even order statistic. The numpy implementation now
+matches other reference implementations.

numpy/lib/_function_base_impl.py

@@ -4355,12 +4355,14 @@ def quantile(a,
     The table above includes only the estimators from H&F that are continuous
     functions of probability `q` (estimators 4-9). NumPy also provides the
     three discontinuous estimators from H&F (estimators 1-3), where ``j`` is
-    defined as above and ``m`` and ``g`` are defined as follows.
+    defined as above, ``m`` is defined as follows, and ``g`` is a function
+    of the real-valued ``index = q*n + m - 1`` and ``j``.
 
-    1. ``inverted_cdf``: ``m = 0`` and ``g = int(q*n > 0)``
-    2. ``averaged_inverted_cdf``: ``m = 0`` and ``g = (1 + int(q*n > 0)) / 2``
+    1. ``inverted_cdf``: ``m = 0`` and ``g = int(index - j > 0)``
+    2. ``averaged_inverted_cdf``: ``m = 0`` and
+       ``g = (1 + int(index - j > 0)) / 2``
     3. ``closest_observation``: ``m = -1/2`` and
-       ``1 - int((g == 0) & (j%2 == 0))``
+       ``g = 1 - int((index == j) & (j%2 == 1))``
 
     For backward compatibility with previous versions of NumPy, `quantile`
     provides four additional discontinuous estimators. Like
@@ -4608,7 +4610,9 @@ def _discret_interpolation_to_boundaries(index, gamma_condition_fun):
 
 
 def _closest_observation(n, quantiles):
-    gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 0)
+    # "choose the nearest even order statistic at g=0" (H&F (1996) pp. 362).
+    # Order is 1-based so for zero-based indexing round to nearest odd index.
+    gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 1)
     return _discret_interpolation_to_boundaries((n * quantiles) - 1 - 0.5,
                                                 gamma_fun)
 

numpy/lib/tests/test_function_base.py

@@ -4029,6 +4029,20 @@ def test_weibull_fraction(self):
         quantile = np.quantile(arr, [Fraction(1, 2)], method='weibull')
         assert_equal(quantile, np.array(Fraction(1, 20)))
 
+    def test_closest_observation(self):
+        # Round ties to nearest even order statistic (see #26656)
+        m = 'closest_observation'
+        q = 0.5
+        arr = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+        assert_equal(2, np.quantile(arr[0:3], q, method=m))
+        assert_equal(2, np.quantile(arr[0:4], q, method=m))
+        assert_equal(2, np.quantile(arr[0:5], q, method=m))
+        assert_equal(3, np.quantile(arr[0:6], q, method=m))
+        assert_equal(4, np.quantile(arr[0:7], q, method=m))
+        assert_equal(4, np.quantile(arr[0:8], q, method=m))
+        assert_equal(4, np.quantile(arr[0:9], q, method=m))
+        assert_equal(5, np.quantile(arr, q, method=m))
+
 
 class TestLerp:
     @hypothesis.given(t0=st.floats(allow_nan=False, allow_infinity=False,