BUG: quantile discrete methods ended up using -1 as index sometimes

Also, the closest-observation did not correctly support multiple
quantiles calculated at the same time (broadcasting error).

Author: seberg

numpy/lib/function_base.py

@@ -4487,19 +4487,22 @@ def _lerp(a, b, t, out=None):
 
 def _get_gamma_mask(shape, default_value, conditioned_value, where):
     out = np.full(shape, default_value)
-    out[where] = conditioned_value
+    np.copyto(out, conditioned_value, where=where, casting="unsafe")
     return out
 
 
 def _discret_interpolation_to_boundaries(index, gamma_condition_fun):
     previous = np.floor(index)
     next = previous + 1
     gamma = index - previous
-    return _get_gamma_mask(shape=index.shape,
-                           default_value=next,
-                           conditioned_value=previous,
-                           where=gamma_condition_fun(gamma, index)
-                           ).astype(np.intp)
+    res = _get_gamma_mask(shape=index.shape,
+                          default_value=next,
+                          conditioned_value=previous,
+                          where=gamma_condition_fun(gamma, index)
+                          ).astype(np.intp)
+    # Some methods can lead to out-of-bound integers, clip them:
+    res[res < 0] = 0
+    return res
 
 
 def _closest_observation(n, quantiles):

numpy/lib/tests/test_function_base.py

@@ -3426,12 +3426,22 @@ def test_quantile_preserve_int_type(self, dtype):
                           method="nearest")
         assert res.dtype == dtype
 
-    def test_quantile_monotonic(self):
+    @pytest.mark.parametrize("method",
+             ['inverted_cdf', 'averaged_inverted_cdf', 'closest_observation',
+              'interpolated_inverted_cdf', 'hazen', 'weibull', 'linear',
+              'median_unbiased', 'normal_unbiased',
+              'nearest', 'lower', 'higher', 'midpoint'])
+    def test_quantile_monotonic(self, method):
         # GH 14685
         # test that the return value of quantile is monotonic if p0 is ordered
-        p0 = np.arange(0, 1, 0.01)
+        # Also tests that the boundary values are not mishandled.
+        p0 = np.linspace(0, 1, 101)
         quantile = np.quantile(np.array([0, 1, 1, 2, 2, 3, 3, 4, 5, 5, 1, 1, 9, 9, 9,
-                                         8, 8, 7]) * 0.1, p0)
+                                         8, 8, 7]) * 0.1, p0, method=method)
+        assert_equal(np.sort(quantile), quantile)
+
+        # Also test one where the number of data points is clearly divisible:
+        quantile = np.quantile([0., 1., 2., 3.], p0, method=method)
         assert_equal(np.sort(quantile), quantile)
 
     @hypothesis.given(