Repository: numpy/numpy
Pull request: BUG,DEP: Fixup quantile/percentile and rename interpolation->method
This (unfortunately currently) does 3 things:

* It renames and deprecates the `interpolation` kwarg as `method`
* It fixes the documentation current plot to use the new methods (see below)
* It fixes two small issues with some discrete methods and tests them by expanding the "monoticity" test case

The main `method` docstring is modified to:
```
    method : str, optional
        This parameter specifies the method to use for estimating the
        percentile.  There are many different methods, some unique to NumPy.
        See the notes for explanation.  The options aligning with the R types
        and the H&F paper are:

        * (H&F 1): 'inverted_cdf'
        * (H&F 2): 'averaged_inverted_cdf'
        * (H&F 3): 'closest_observation'
        * (H&F 4): 'interpolated_inverted_cdf'
        * (H&F 5): 'hazen'
        * (H&F 6): 'weibull'
        * (H&F 7): 'linear'  (default)
        * (H&F 8): 'median_unbiased'
        * (H&F 9): 'normal_unbiased'

        Mainly for compatibility reasons, NumPy also supports the following
        options which appear to be unique to NumPy:

        * 'lower'
        * 'higher',
        * 'midpoint'
        * 'nearest'

        .. versionchanged:: 1.22.0
            This argument was previously called "interpolation" and only
            offered the "linear" default and last four options.
```

Here is the plot, note that it is smaller in the artifact: https://23201-908607-gh.circle-artifacts.com/0/doc/build/html/reference/generated/numpy.percentile.html

![new-percentile-plot](https://user-images.githubusercontent.com/61977/140831414-b88a4954-fcd2-42fd-a79d-f92422d8b665.png)
(yeah, colors/linestyles are a bit annoying, some of the colors match, because they are somewhat related, I think)

Closes gh-20306.

Problem statement:
DEP: Rename quantile/percentile `interpolation` kwarg to `method`
It seems like a good idea to rename the quantile/percentile "interpolation" kwarg to method.  This is because:
* Method seems more descriptive, "interpolation" doesn't actually capture all of the complexity (anymore)
* The old versions (except the default) were not really described in literature anyway.  So this nudges users to switch to a different value.

Hints from issue discussion:

Commit message: DOC: Add release not for quantile `interpolation` rename

Also updates the old release note to include the info (and generally
tweaks it a bit)

Patch to review:
diff --git a/doc/release/upcoming_changes/19857.improvement.rst b/doc/release/upcoming_changes/19857.improvement.rst
index e39d413cc3bc..cbeff08b1ade 100644
--- a/doc/release/upcoming_changes/19857.improvement.rst
+++ b/doc/release/upcoming_changes/19857.improvement.rst
@@ -1,8 +1,13 @@
-Add new linear interpolation methods for ``quantile`` and ``percentile``
-------------------------------------------------------------------------
+Add new methods for ``quantile`` and ``percentile``
+---------------------------------------------------
 
-``quantile`` and ``percentile`` now have 13 linear interpolation methods,
-nine of which can be found in the scientific literature.
-The remaining methods are NumPy specific and are kept for backwards
-compatibility. The default is "inclusive" (method 7), whose behavior is equivalent
-to the former default "linear".
+``quantile`` and ``percentile`` now have have a ``method=``
+keyword argument supporting 13 different methods.
+This replaces the ``interpolation=`` keyword argument.
+
+The methods are now aligned with nine methods which can be
+found in scientific literature and the R language.
+The remaining methods are the previous discontinuous variations
+of the default "linear" one.
+
+Please see the documentation of `numpy.percentile` for more information.
diff --git a/numpy/core/tests/test_deprecations.py b/numpy/core/tests/test_deprecations.py
index a1b379d9274d..93e07216b2f8 100644
--- a/numpy/core/tests/test_deprecations.py
+++ b/numpy/core/tests/test_deprecations.py
@@ -1230,3 +1230,23 @@ def test_deprecated_module(self):
     def test_deprecated_attr(self):
         finfo = np.finfo(float)
         self.assert_deprecated(lambda: getattr(finfo, "machar"))
+
+
+class TestQuantileInterpolationDeprecation(_DeprecationTestCase):
+    # Deprecated 2021-11-08, NumPy 1.22
+    @pytest.mark.parametrize("func",
+        [np.percentile, np.quantile, np.nanpercentile, np.nanquantile])
+    def test_deprecated(self, func):
+        self.assert_deprecated(
+            lambda: func([0., 1.], 0., interpolation="linear"))
+        self.assert_deprecated(
+            lambda: func([0., 1.], 0., interpolation="nearest"))
+
+    @pytest.mark.parametrize("func",
+            [np.percentile, np.quantile, np.nanpercentile, np.nanquantile])
+    def test_both_passed(self, func):
+        with warnings.catch_warnings():
+            # catch the warning, but make sure it does not raise:
+            warnings.simplefilter("always", DeprecationWarning)
+            with pytest.raises(TypeError):
+                func([0., 1.], 0., interpolation="nearest", method="nearest")
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index 3c9983edf31f..6d84627cdb09 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -50,8 +50,8 @@
     'quantile'
     ]
 
-# _QuantileInterpolation is a dictionary listing all the supported
-# interpolation methods to compute quantile/percentile.
+# _QuantileMethods is a dictionary listing all the supported methods to
+# compute quantile/percentile.
 #
 # Below virtual_index refer to the index of the element where the percentile
 # would be found in the sorted sample.
@@ -61,13 +61,13 @@
 # is made of a integer part (a.k.a 'i' or 'left') and a fractional part
 # (a.k.a 'g' or 'gamma')
 #
-# Each _QuantileInterpolation have two properties
+# Each _QuantileMethods has two properties
 # get_virtual_index : Callable
 #   The function used to compute the virtual_index.
 # fix_gamma : Callable
 #   A function used for discret methods to force the index to a specific value.
-_QuantileInterpolation = dict(
-    # --- HYNDMAN AND FAN METHODS
+_QuantileMethods = dict(
+    # --- HYNDMAN and FAN METHODS
     # Discrete methods
     inverted_cdf=dict(
         get_virtual_index=lambda n, quantiles: _inverted_cdf(n, quantiles),
@@ -3854,7 +3854,7 @@ def _median(a, axis=None, out=None, overwrite_input=False):
 
 
 def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
-                           interpolation=None, keepdims=None):
+                           method=None, keepdims=None, *, interpolation=None):
     return (a, q, out)
 
 
@@ -3864,8 +3864,10 @@ def percentile(a,
                axis=None,
                out=None,
                overwrite_input=False,
-               interpolation="linear",
-               keepdims=False):
+               method="linear",
+               keepdims=False,
+               *,
+               interpolation=None):
     """
     Compute the q-th percentile of the data along the specified axis.
 
@@ -3893,31 +3895,33 @@ def percentile(a,
         If True, then allow the input array `a` to be modified by intermediate
         calculations, to save memory. In this case, the contents of the input
         `a` after this function completes is undefined.
-    interpolation : str, optional
-        This parameter specifies the interpolation method to
-        use when the desired percentile lies between two data points
-        There are many different methods, some unique to NumPy. See the
-        notes for explanation. Options
-
-        * (NPY 1): 'lower'
-        * (NPY 2): 'higher',
-        * (NPY 3): 'midpoint'
-        * (NPY 4): 'nearest'
-        * (NPY 5): 'linear'
-
-        New options:
-
-        * (H&F 1): 'inverted_cdf'
-        * (H&F 2): 'averaged_inverted_cdf'
-        * (H&F 3): 'closest_observation'
-        * (H&F 4): 'interpolated_inverted_cdf'
-        * (H&F 5): 'hazen'
-        * (H&F 6): 'weibull'
-        * (H&F 7): 'linear'  (default)
-        * (H&F 8): 'median_unbiased'
-        * (H&F 9): 'normal_unbiased'
+    method : str, optional
+        This parameter specifies the method to use for estimating the
+        percentile.  There are many different methods, some unique to NumPy.
+        See the notes for explanation.  The options sorted by their R type
+        as summarized in the H&F paper [1]_ are:
+
+        1. 'inverted_cdf'
+        2. 'averaged_inverted_cdf'
+        3. 'closest_observation'
+        4. 'interpolated_inverted_cdf'
+        5. 'hazen'
+        6. 'weibull'
+        7. 'linear'  (default)
+        8. 'median_unbiased'
+        9. 'normal_unbiased'
+
+        The first three methods are discontiuous.  NumPy further defines the
+        following discontinuous variations of the default 'linear' (7.) option:
+
+        * 'lower'
+        * 'higher',
+        * 'midpoint'
+        * 'nearest'
 
         .. versionchanged:: 1.22.0
+            This argument was previously called "interpolation" and only
+            offered the "linear" default and last four options.
 
     keepdims : bool, optional
         If this is set to True, the axes which are reduced are left in
@@ -3926,6 +3930,11 @@ def percentile(a,
 
         .. versionadded:: 1.9.0
 
+    interpolation : str, optional
+        Deprecated name for the method keyword argument.
+
+        .. deprecated:: 1.22.0
+
     Returns
     -------
     percentile : scalar or ndarray
@@ -3950,16 +3959,16 @@ def percentile(a,
     Given a vector ``V`` of length ``N``, the q-th percentile of ``V`` is
     the value ``q/100`` of the way from the minimum to the maximum in a
     sorted copy of ``V``. The values and distances of the two nearest
-    neighbors as well as the `interpolation` parameter will determine the
+    neighbors as well as the `method` parameter will determine the
     percentile if the normalized ranking does not match the location of
     ``q`` exactly. This function is the same as the median if ``q=50``, the
     same as the minimum if ``q=0`` and the same as the maximum if
     ``q=100``.
 
-    This optional `interpolation` parameter specifies the interpolation
-    method to use when the desired quantile lies between two data points
-    ``i < j``. If ``g`` is the fractional part of the index surrounded by
-    ``i`` and  alpha and beta are correction constants modifying i and j.
+    This optional `method` parameter specifies the method to use when the
+    desired quantile lies between two data points ``i < j``.
+    If ``g`` is the fractional part of the index surrounded by ``i`` and
+    alpha and beta are correction constants modifying i and j.
 
     Below, 'q' is the quantile value, 'n' is the sample size and
     alpha and beta are constants.
@@ -3970,7 +3979,7 @@ def percentile(a,
     .. math::
         i + g = (q - alpha) / ( n - alpha - beta + 1 )
 
-    The different interpolation methods then work as follows
+    The different methods then work as follows
 
     inverted_cdf:
         method 1 of H&F [1]_.
@@ -4075,7 +4084,7 @@ def percentile(a,
     array([7.,  2.])
     >>> assert not np.all(a == b)
 
-    The different types of interpolation can be visualized graphically:
+    The different methods can be visualized graphically:
 
     .. plot::
 
@@ -4085,20 +4094,25 @@ def percentile(a,
         p = np.linspace(0, 100, 6001)
         ax = plt.gca()
         lines = [
-            ('linear', None),
-            ('higher', '--'),
-            ('lower', '--'),
-            ('nearest', '-.'),
-            ('midpoint', '-.'),
-        ]
-        for interpolation, style in lines:
+            ('linear', '-', 'C0'),
+            ('inverted_cdf', ':', 'C1'),
+            # Almost the same as `inverted_cdf`:
+            ('averaged_inverted_cdf', '-.', 'C1'),
+            ('closest_observation', ':', 'C2'),
+            ('interpolated_inverted_cdf', '--', 'C1'),
+            ('hazen', '--', 'C3'),
+            ('weibull', '-.', 'C4'),
+            ('median_unbiased', '--', 'C5'),
+            ('normal_unbiased', '-.', 'C6'),
+            ]
+        for method, style, color in lines:
             ax.plot(
-                p, np.percentile(a, p, interpolation=interpolation),
-                label=interpolation, linestyle=style)
+                p, np.percentile(a, p, method=method),
+                label=method, linestyle=style, color=color)
         ax.set(
-            title='Interpolation methods for list: ' + str(a),
+            title='Percentiles for different methods and data: ' + str(a),
             xlabel='Percentile',
-            ylabel='List item returned',
+            ylabel='Estimated percentile value',
             yticks=a)
         ax.legend()
         plt.show()
@@ -4110,16 +4124,19 @@ def percentile(a,
        The American Statistician, 50(4), pp. 361-365, 1996
 
     """
+    if interpolation is not None:
+        method = _check_interpolation_as_method(
+            method, interpolation, "percentile")
     q = np.true_divide(q, 100)
     q = asanyarray(q)  # undo any decay that the ufunc performed (see gh-13105)
     if not _quantile_is_valid(q):
         raise ValueError("Percentiles must be in the range [0, 100]")
     return _quantile_unchecked(
-        a, q, axis, out, overwrite_input, interpolation, keepdims)
+        a, q, axis, out, overwrite_input, method, keepdims)
 
 
 def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
-                         interpolation=None, keepdims=None):
+                         method=None, keepdims=None, *, interpolation=None):
     return (a, q, out)
 
 
@@ -4129,8 +4146,10 @@ def quantile(a,
              axis=None,
              out=None,
              overwrite_input=False,
-             interpolation="linear",
-             keepdims=False):
+             method="linear",
+             keepdims=False,
+             *,
+             interpolation=None):
     """
     Compute the q-th quantile of the data along the specified axis.
 
@@ -4155,37 +4174,44 @@ def quantile(a,
         intermediate calculations, to save memory. In this case, the
         contents of the input `a` after this function completes is
         undefined.
-    interpolation : str, optional
-        This parameter specifies the interpolation method to use when the
-        desired quantile lies between two data points There are many
-        different methods, some unique to NumPy. See the notes for
-        explanation. Options:
-
-        * (NPY 1): 'lower'
-        * (NPY 2): 'higher',
-        * (NPY 3): 'midpoint'
-        * (NPY 4): 'nearest'
-        * (NPY 5): 'linear'
-
-        New options:
-
-        * (H&F 1): 'inverted_cdf'
-        * (H&F 2): 'averaged_inverted_cdf'
-        * (H&F 3): 'closest_observation'
-        * (H&F 4): 'interpolated_inverted_cdf'
-        * (H&F 5): 'hazen'
-        * (H&F 6): 'weibull'
-        * (H&F 7): 'linear'  (default)
-        * (H&F 8): 'median_unbiased'
-        * (H&F 9): 'normal_unbiased'
-
-        .. versionadded:: 1.22.0
+    method : str, optional
+        This parameter specifies the method to use for estimating the
+        quantile.  There are many different methods, some unique to NumPy.
+        See the notes for explanation.  The options sorted by their R type
+        as summarized in the H&F paper [1]_ are:
+
+        1. 'inverted_cdf'
+        2. 'averaged_inverted_cdf'
+        3. 'closest_observation'
+        4. 'interpolated_inverted_cdf'
+        5. 'hazen'
+        6. 'weibull'
+        7. 'linear'  (default)
+        8. 'median_unbiased'
+        9. 'normal_unbiased'
+
+        The first three methods are discontiuous.  NumPy further defines the
+        following discontinuous variations of the default 'linear' (7.) option:
+
+        * 'lower'
+        * 'higher',
+        * 'midpoint'
+        * 'nearest'
+
+        .. versionchanged:: 1.22.0
+            This argument was previously called "interpolation" and only
+            offered the "linear" default and last four options.
 
     keepdims : bool, optional
         If this is set to True, the axes which are reduced are left in
         the result as dimensions with size one. With this option, the
         result will broadcast correctly against the original array `a`.
 
+    interpolation : str, optional
+        Deprecated name for the method keyword argument.
+
+        .. deprecated:: 1.22.0
+
     Returns
     -------
     quantile : scalar or ndarray
@@ -4210,20 +4236,20 @@ def quantile(a,
     Given a vector ``V`` of length ``N``, the q-th quantile of ``V`` is the
     value ``q`` of the way from the minimum to the maximum in a sorted copy of
     ``V``. The values and distances of the two nearest neighbors as well as the
-    `interpolation` parameter will determine the quantile if the normalized
+    `method` parameter will determine the quantile if the normalized
     ranking does not match the location of ``q`` exactly. This function is the
     same as the median if ``q=0.5``, the same as the minimum if ``q=0.0`` and
     the same as the maximum if ``q=1.0``.
 
-    This optional `interpolation` parameter specifies the interpolation method
-    to use when the desired quantile lies between two data points ``i < j``. If
-    ``g`` is the fractional part of the index surrounded by ``i`` and  alpha
-    and beta are correction constants modifying i and j.
+    This optional `method` parameter specifies the method to use when the
+    desired quantile lies between two data points ``i < j``.
+    If ``g`` is the fractional part of the index surrounded by ``i`` and
+    alpha and beta are correction constants modifying i and j.
 
     .. math::
         i + g = (q - alpha) / ( n - alpha - beta + 1 )
 
-    The different interpolation methods then work as follows
+    The different methods then work as follows
 
     inverted_cdf:
         method 1 of H&F [1]_.
@@ -4326,6 +4352,8 @@ def quantile(a,
     array([7.,  2.])
     >>> assert not np.all(a == b)
 
+    See also `numpy.percentile` for a visualization of most methods.
+
     References
     ----------
     .. [1] R. J. Hyndman and Y. Fan,
@@ -4333,11 +4361,15 @@ def quantile(a,
        The American Statistician, 50(4), pp. 361-365, 1996
 
     """
+    if interpolation is not None:
+        method = _check_interpolation_as_method(
+            method, interpolation, "quantile")
+
     q = np.asanyarray(q)
     if not _quantile_is_valid(q):
         raise ValueError("Quantiles must be in the range [0, 1]")
     return _quantile_unchecked(
-        a, q, axis, out, overwrite_input, interpolation, keepdims)
+        a, q, axis, out, overwrite_input, method, keepdims)
 
 
 def _quantile_unchecked(a,
@@ -4345,7 +4377,7 @@ def _quantile_unchecked(a,
                         axis=None,
                         out=None,
                         overwrite_input=False,
-                        interpolation="linear",
+                        method="linear",
                         keepdims=False):
     """Assumes that q is in [0, 1], and is an ndarray"""
     r, k = _ureduce(a,
@@ -4354,7 +4386,7 @@ def _quantile_unchecked(a,
                     axis=axis,
                     out=out,
                     overwrite_input=overwrite_input,
-                    interpolation=interpolation)
+                    method=method)
     if keepdims:
         return r.reshape(q.shape + k)
     else:
@@ -4373,6 +4405,23 @@ def _quantile_is_valid(q):
     return True
 
 
+def _check_interpolation_as_method(method, interpolation, fname):
+    # Deprecated NumPy 1.22, 2021-11-08
+    warnings.warn(
+        f"the `interpolation=` argument to {fname} was renamed to "
+        "`method=`, which has additional options.\n"
+        "Users of the modes 'nearest', 'lower', 'higher', or "
+        "'midpoint' are encouraged to review the method they. "
+        "(Deprecated NumPy 1.22)",
+        DeprecationWarning, stacklevel=4)
+    if method != "linear":
+        # sanity check, we assume this basically never happens
+        raise TypeError(
+            "You shall not pass both `method` and `interpolation`!\n"
+            "(`interpolation` is Deprecated in favor of `method`)")
+    return interpolation
+
+
 def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
     """
     Compute the floating point indexes of an array for the linear
@@ -4398,9 +4447,7 @@ def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
     ) - 1
 
 
-def _get_gamma(virtual_indexes,
-               previous_indexes,
-               interpolation: _QuantileInterpolation):
+def _get_gamma(virtual_indexes, previous_indexes, method):
     """
     Compute gamma (a.k.a 'm' or 'weight') for the linear interpolation
     of quantiles.
@@ -4410,7 +4457,7 @@ def _get_gamma(virtual_indexes,
         sample.
     previous_indexes : array_like
         The floor values of virtual_indexes.
-    interpolation : _QuantileInterpolation
+    interpolation : dict
         The interpolation method chosen, which may have a specific rule
         modifying gamma.
 
@@ -4418,7 +4465,7 @@ def _get_gamma(virtual_indexes,
     by the interpolation method.
     """
     gamma = np.asanyarray(virtual_indexes - previous_indexes)
-    gamma = interpolation["fix_gamma"](gamma, virtual_indexes)
+    gamma = method["fix_gamma"](gamma, virtual_indexes)
     return np.asanyarray(gamma)
 
 
@@ -4447,7 +4494,7 @@ 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
 
 
@@ -4455,11 +4502,14 @@ 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):
@@ -4480,7 +4530,7 @@ def _quantile_ureduce_func(
         axis: int = None,
         out=None,
         overwrite_input: bool = False,
-        interpolation="linear",
+        method="linear",
 ) -> np.array:
     if q.ndim > 2:
         # The code below works fine for nd, but it might not have useful
@@ -4502,7 +4552,7 @@ def _quantile_ureduce_func(
     result = _quantile(arr,
                        quantiles=q,
                        axis=axis,
-                       interpolation=interpolation,
+                       method=method,
                        out=out)
     return result
 
@@ -4546,7 +4596,7 @@ def _quantile(
         arr: np.array,
         quantiles: np.array,
         axis: int = -1,
-        interpolation="linear",
+        method="linear",
         out=None,
 ):
     """
@@ -4556,8 +4606,8 @@ def _quantile(
     It computes the quantiles of the array for the given axis.
     A linear interpolation is performed based on the `interpolation`.
 
-    By default, the interpolation is "linear" where
-    alpha == beta == 1 which performs the 7th method of Hyndman&Fan.
+    By default, the method is "linear" where alpha == beta == 1 which
+    performs the 7th method of Hyndman&Fan.
     With "median_unbiased" we get alpha == beta == 1/3
     thus the 8th method of Hyndman&Fan.
     """
@@ -4574,13 +4624,12 @@ def _quantile(
     # Virtual because it is a floating point value, not an valid index.
     # The nearest neighbours are used for interpolation
     try:
-        interpolation = _QuantileInterpolation[interpolation]
+        method = _QuantileMethods[method]
     except KeyError:
         raise ValueError(
-            f"{interpolation!r} is not a valid interpolation. Use one of: "
-            f"{_QuantileInterpolation.keys()}") from None
-    virtual_indexes = interpolation["get_virtual_index"](values_count,
-                                                         quantiles)
+            f"{method!r} is not a valid method. Use one of: "
+            f"{_QuantileMethods.keys()}") from None
+    virtual_indexes = method["get_virtual_index"](values_count, quantiles)
     virtual_indexes = np.asanyarray(virtual_indexes)
     if np.issubdtype(virtual_indexes.dtype, np.integer):
         # No interpolation needed, take the points along axis
@@ -4614,9 +4663,7 @@ def _quantile(
         previous = np.take(arr, previous_indexes, axis=DATA_AXIS)
         next = np.take(arr, next_indexes, axis=DATA_AXIS)
         # --- Linear interpolation
-        gamma = _get_gamma(virtual_indexes,
-                           previous_indexes,
-                           interpolation)
+        gamma = _get_gamma(virtual_indexes, previous_indexes, method)
         result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
         gamma = gamma.reshape(result_shape)
         result = _lerp(previous,
diff --git a/numpy/lib/function_base.pyi b/numpy/lib/function_base.pyi
index 82c625fed9d6..7e227f9da52d 100644
--- a/numpy/lib/function_base.pyi
+++ b/numpy/lib/function_base.pyi
@@ -500,7 +500,7 @@ def median(
     keepdims: bool = ...,
 ) -> _ArrayType: ...
 
-_InterpolationKind = L[
+_MethodKind = L[
     "inverted_cdf",
     "averaged_inverted_cdf",
     "closest_observation",
@@ -523,7 +523,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> floating[Any]: ...
 @overload
@@ -533,7 +533,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> complexfloating[Any, Any]: ...
 @overload
@@ -543,7 +543,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> timedelta64: ...
 @overload
@@ -553,7 +553,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> datetime64: ...
 @overload
@@ -563,7 +563,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> Any: ...
 @overload
@@ -573,7 +573,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> NDArray[floating[Any]]: ...
 @overload
@@ -583,7 +583,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> NDArray[complexfloating[Any, Any]]: ...
 @overload
@@ -593,7 +593,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> NDArray[timedelta64]: ...
 @overload
@@ -603,7 +603,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> NDArray[datetime64]: ...
 @overload
@@ -613,7 +613,7 @@ def percentile(
     axis: None = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: L[False] = ...,
 ) -> NDArray[object_]: ...
 @overload
@@ -623,7 +623,7 @@ def percentile(
     axis: None | _ShapeLike = ...,
     out: None = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: bool = ...,
 ) -> Any: ...
 @overload
@@ -633,7 +633,7 @@ def percentile(
     axis: None | _ShapeLike = ...,
     out: _ArrayType = ...,
     overwrite_input: bool = ...,
-    interpolation: _InterpolationKind = ...,
+    method: _MethodKind = ...,
     keepdims: bool = ...,
 ) -> _ArrayType: ...
 
diff --git a/numpy/lib/nanfunctions.py b/numpy/lib/nanfunctions.py
index 7e953be03757..d7ea1ca65b7d 100644
--- a/numpy/lib/nanfunctions.py
+++ b/numpy/lib/nanfunctions.py
@@ -1223,8 +1223,9 @@ def nanmedian(a, axis=None, out=None, overwrite_input=False, keepdims=np._NoValu
         return r
 
 
-def _nanpercentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
-                              interpolation=None, keepdims=None):
+def _nanpercentile_dispatcher(
+        a, q, axis=None, out=None, overwrite_input=None,
+        method=None, keepdims=None, *, interpolation=None):
     return (a, q, out)
 
 
@@ -1235,8 +1236,10 @@ def nanpercentile(
         axis=None,
         out=None,
         overwrite_input=False,
-        interpolation="linear",
+        method="linear",
         keepdims=np._NoValue,
+        *,
+        interpolation=None,
 ):
     """
     Compute the qth percentile of the data along the specified axis,
@@ -1267,31 +1270,33 @@ def nanpercentile(
         intermediate calculations, to save memory. In this case, the
         contents of the input `a` after this function completes is
         undefined.
-    interpolation : str, optional
-        This parameter specifies the interpolation method to use when the
-        desired percentile lies between two data points There are many
-        different methods, some unique to NumPy. See the notes for
-        explanation. Options:
-
-        * (NPY 1): 'lower'
-        * (NPY 2): 'higher',
-        * (NPY 3): 'midpoint'
-        * (NPY 4): 'nearest'
-        * (NPY 5): 'linear'  (default)
-
-        New options:
-
-        * (H&F 1): 'inverted_cdf'
-        * (H&F 2): 'averaged_inverted_cdf'
-        * (H&F 3): 'closest_observation'
-        * (H&F 4): 'interpolated_inverted_cdf'
-        * (H&F 5): 'hazen'
-        * (H&F 6): 'weibull'
-        * (H&F 7): 'linear'  (default)
-        * (H&F 8): 'median_unbiased'
-        * (H&F 9): 'normal_unbiased'
+    method : str, optional
+        This parameter specifies the method to use for estimating the
+        percentile.  There are many different methods, some unique to NumPy.
+        See the notes for explanation.  The options sorted by their R type
+        as summarized in the H&F paper [1]_ are:
+
+        1. 'inverted_cdf'
+        2. 'averaged_inverted_cdf'
+        3. 'closest_observation'
+        4. 'interpolated_inverted_cdf'
+        5. 'hazen'
+        6. 'weibull'
+        7. 'linear'  (default)
+        8. 'median_unbiased'
+        9. 'normal_unbiased'
+
+        The first three methods are discontiuous.  NumPy further defines the
+        following discontinuous variations of the default 'linear' (7.) option:
+
+        * 'lower'
+        * 'higher',
+        * 'midpoint'
+        * 'nearest'
 
-        .. versionadded:: 1.22.0
+        .. versionchanged:: 1.22.0
+            This argument was previously called "interpolation" and only
+            offered the "linear" default and last four options.
 
     keepdims : bool, optional
         If this is set to True, the axes which are reduced are left in
@@ -1304,6 +1309,11 @@ def nanpercentile(
         a sub-class and `mean` does not have the kwarg `keepdims` this
         will raise a RuntimeError.
 
+    interpolation : str, optional
+        Deprecated name for the method keyword argument.
+
+        .. deprecated:: 1.22.0
+
     Returns
     -------
     percentile : scalar or ndarray
@@ -1355,7 +1365,17 @@ def nanpercentile(
     array([7., 2.])
     >>> assert not np.all(a==b)
 
+    References
+    ----------
+    .. [1] R. J. Hyndman and Y. Fan,
+       "Sample quantiles in statistical packages,"
+       The American Statistician, 50(4), pp. 361-365, 1996
+
     """
+    if interpolation is not None:
+        method = function_base._check_interpolation_as_method(
+            method, interpolation, "nanpercentile")
+
     a = np.asanyarray(a)
     q = np.true_divide(q, 100.0)
     # undo any decay that the ufunc performed (see gh-13105)
@@ -1363,11 +1383,11 @@ def nanpercentile(
     if not function_base._quantile_is_valid(q):
         raise ValueError("Percentiles must be in the range [0, 100]")
     return _nanquantile_unchecked(
-        a, q, axis, out, overwrite_input, interpolation, keepdims)
+        a, q, axis, out, overwrite_input, method, keepdims)
 
 
 def _nanquantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
-                            interpolation=None, keepdims=None):
+                            method=None, keepdims=None, *, interpolation=None):
     return (a, q, out)
 
 
@@ -1378,8 +1398,10 @@ def nanquantile(
         axis=None,
         out=None,
         overwrite_input=False,
-        interpolation="linear",
+        method="linear",
         keepdims=np._NoValue,
+        *,
+        interpolation=None,
 ):
     """
     Compute the qth quantile of the data along the specified axis,
@@ -1408,31 +1430,33 @@ def nanquantile(
         If True, then allow the input array `a` to be modified by intermediate
         calculations, to save memory. In this case, the contents of the input
         `a` after this function completes is undefined.
-    interpolation : str, optional
-        This parameter specifies the interpolation method to
-        use when the desired quantile lies between two data points
-        There are many different methods, some unique to NumPy. See the
-        notes for explanation. Options:
-
-        * (NPY 1): 'lower'
-        * (NPY 2): 'higher',
-        * (NPY 3): 'midpoint'
-        * (NPY 4): 'nearest'
-        * (NPY 5): 'linear'  (default)
-
-        New options:
-
-        * (H&F 1): 'inverted_cdf'
-        * (H&F 2): 'averaged_inverted_cdf'
-        * (H&F 3): 'closest_observation'
-        * (H&F 4): 'interpolated_inverted_cdf'
-        * (H&F 5): 'hazen'
-        * (H&F 6): 'weibull'
-        * (H&F 7): 'linear'  (default)
-        * (H&F 8): 'median_unbiased'
-        * (H&F 9): 'normal_unbiased'
+    method : str, optional
+        This parameter specifies the method to use for estimating the
+        quantile.  There are many different methods, some unique to NumPy.
+        See the notes for explanation.  The options sorted by their R type
+        as summarized in the H&F paper [1]_ are:
+
+        1. 'inverted_cdf'
+        2. 'averaged_inverted_cdf'
+        3. 'closest_observation'
+        4. 'interpolated_inverted_cdf'
+        5. 'hazen'
+        6. 'weibull'
+        7. 'linear'  (default)
+        8. 'median_unbiased'
+        9. 'normal_unbiased'
+
+        The first three methods are discontiuous.  NumPy further defines the
+        following discontinuous variations of the default 'linear' (7.) option:
+
+        * 'lower'
+        * 'higher',
+        * 'midpoint'
+        * 'nearest'
 
         .. versionchanged:: 1.22.0
+            This argument was previously called "interpolation" and only
+            offered the "linear" default and last four options.
 
     keepdims : bool, optional
         If this is set to True, the axes which are reduced are left in
@@ -1445,6 +1469,11 @@ def nanquantile(
         a sub-class and `mean` does not have the kwarg `keepdims` this
         will raise a RuntimeError.
 
+    interpolation : str, optional
+        Deprecated name for the method keyword argument.
+
+        .. deprecated:: 1.22.0
+
     Returns
     -------
     quantile : scalar or ndarray
@@ -1495,13 +1524,23 @@ def nanquantile(
     array([7., 2.])
     >>> assert not np.all(a==b)
 
+    References
+    ----------
+    .. [1] R. J. Hyndman and Y. Fan,
+       "Sample quantiles in statistical packages,"
+       The American Statistician, 50(4), pp. 361-365, 1996
+
     """
+    if interpolation is not None:
+        method = function_base._check_interpolation_as_method(
+            method, interpolation, "nanquantile")
+
     a = np.asanyarray(a)
     q = np.asanyarray(q)
     if not function_base._quantile_is_valid(q):
         raise ValueError("Quantiles must be in the range [0, 1]")
     return _nanquantile_unchecked(
-        a, q, axis, out, overwrite_input, interpolation, keepdims)
+        a, q, axis, out, overwrite_input, method, keepdims)
 
 
 def _nanquantile_unchecked(
@@ -1510,7 +1549,7 @@ def _nanquantile_unchecked(
         axis=None,
         out=None,
         overwrite_input=False,
-        interpolation="linear",
+        method="linear",
         keepdims=np._NoValue,
 ):
     """Assumes that q is in [0, 1], and is an ndarray"""
@@ -1524,7 +1563,7 @@ def _nanquantile_unchecked(
                                   axis=axis,
                                   out=out,
                                   overwrite_input=overwrite_input,
-                                  interpolation=interpolation)
+                                  method=method)
     if keepdims and keepdims is not np._NoValue:
         return r.reshape(q.shape + k)
     else:
@@ -1532,7 +1571,7 @@ def _nanquantile_unchecked(
 
 
 def _nanquantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False,
-                              interpolation="linear"):
+                              method="linear"):
     """
     Private function that doesn't support extended axis or keepdims.
     These methods are extended to this function using _ureduce
@@ -1540,10 +1579,10 @@ def _nanquantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False,
     """
     if axis is None or a.ndim == 1:
         part = a.ravel()
-        result = _nanquantile_1d(part, q, overwrite_input, interpolation)
+        result = _nanquantile_1d(part, q, overwrite_input, method)
     else:
         result = np.apply_along_axis(_nanquantile_1d, axis, a, q,
-                                     overwrite_input, interpolation)
+                                     overwrite_input, method)
         # apply_along_axis fills in collapsed axis with results.
         # Move that axis to the beginning to match percentile's
         # convention.
@@ -1555,7 +1594,7 @@ def _nanquantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False,
     return result
 
 
-def _nanquantile_1d(arr1d, q, overwrite_input=False, interpolation="linear"):
+def _nanquantile_1d(arr1d, q, overwrite_input=False, method="linear"):
     """
     Private function for rank 1 arrays. Compute quantile ignoring NaNs.
     See nanpercentile for parameter usage
@@ -1567,7 +1606,7 @@ def _nanquantile_1d(arr1d, q, overwrite_input=False, interpolation="linear"):
         return np.full(q.shape, np.nan, dtype=arr1d.dtype)[()]
 
     return function_base._quantile_unchecked(
-        arr1d, q, overwrite_input=overwrite_input, interpolation=interpolation)
+        arr1d, q, overwrite_input=overwrite_input, method=method)
 
 
 def _nanvar_dispatcher(a, axis=None, dtype=None, out=None, ddof=None,
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index 1c274afae2f7..b67a31b1850e 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -2867,7 +2867,7 @@ def test_basic(self):
         assert_equal(np.percentile(x, 50), 1.75)
         x[1] = np.nan
         assert_equal(np.percentile(x, 0), np.nan)
-        assert_equal(np.percentile(x, 0, interpolation='nearest'), np.nan)
+        assert_equal(np.percentile(x, 0, method='nearest'), np.nan)
 
     def test_fraction(self):
         x = [Fraction(i, 2) for i in range(8)]
@@ -2910,7 +2910,7 @@ def test_linear_nan_1D(self, dtype):
         res = np.percentile(
             arr,
             40.0,
-            interpolation="linear")
+            method="linear")
         np.testing.assert_equal(res, np.NAN)
         np.testing.assert_equal(res.dtype, arr.dtype)
 
@@ -2926,7 +2926,7 @@ def test_linear_nan_1D(self, dtype):
                            (np.dtype("O"), np.float64)]
 
     @pytest.mark.parametrize(["input_dtype", "expected_dtype"], H_F_TYPE_CODES)
-    @pytest.mark.parametrize(["interpolation", "expected"],
+    @pytest.mark.parametrize(["method", "expected"],
                              [("inverted_cdf", 20),
                               ("averaged_inverted_cdf", 27.5),
                               ("closest_observation", 20),
@@ -2938,16 +2938,16 @@ def test_linear_nan_1D(self, dtype):
                               ("normal_unbiased", 27.125),
                               ])
     def test_linear_interpolation(self,
-                                  interpolation,
+                                  method,
                                   expected,
                                   input_dtype,
                                   expected_dtype):
         arr = np.asarray([15.0, 20.0, 35.0, 40.0, 50.0], dtype=input_dtype)
-        actual = np.percentile(arr, 40.0, interpolation=interpolation)
+        actual = np.percentile(arr, 40.0, method=method)
 
         np.testing.assert_almost_equal(actual, expected, 14)
 
-        if interpolation in ["inverted_cdf", "closest_observation"]:
+        if method in ["inverted_cdf", "closest_observation"]:
             if input_dtype == "O":
                 np.testing.assert_equal(np.asarray(actual).dtype, np.float64)
             else:
@@ -2962,27 +2962,27 @@ def test_linear_interpolation(self,
     @pytest.mark.parametrize("dtype", TYPE_CODES)
     def test_lower_higher(self, dtype):
         assert_equal(np.percentile(np.arange(10, dtype=dtype), 50,
-                                   interpolation='lower'), 4)
+                                   method='lower'), 4)
         assert_equal(np.percentile(np.arange(10, dtype=dtype), 50,
-                                   interpolation='higher'), 5)
+                                   method='higher'), 5)
 
     @pytest.mark.parametrize("dtype", TYPE_CODES)
     def test_midpoint(self, dtype):
         assert_equal(np.percentile(np.arange(10, dtype=dtype), 51,
-                                   interpolation='midpoint'), 4.5)
+                                   method='midpoint'), 4.5)
         assert_equal(np.percentile(np.arange(9, dtype=dtype) + 1, 50,
-                                   interpolation='midpoint'), 5)
+                                   method='midpoint'), 5)
         assert_equal(np.percentile(np.arange(11, dtype=dtype), 51,
-                                   interpolation='midpoint'), 5.5)
+                                   method='midpoint'), 5.5)
         assert_equal(np.percentile(np.arange(11, dtype=dtype), 50,
-                                   interpolation='midpoint'), 5)
+                                   method='midpoint'), 5)
 
     @pytest.mark.parametrize("dtype", TYPE_CODES)
     def test_nearest(self, dtype):
         assert_equal(np.percentile(np.arange(10, dtype=dtype), 51,
-                                   interpolation='nearest'), 5)
+                                   method='nearest'), 5)
         assert_equal(np.percentile(np.arange(10, dtype=dtype), 49,
-                                       interpolation='nearest'), 4)
+                                   method='nearest'), 4)
 
     def test_linear_interpolation_extrapolation(self):
         arr = np.random.rand(5)
@@ -3019,19 +3019,19 @@ def test_axis(self):
         assert_equal(
             np.percentile(x, (25, 50, 75), axis=1).shape, (3, 3, 5, 6))
         assert_equal(np.percentile(x, (25, 50),
-                                   interpolation="higher").shape, (2,))
+                                   method="higher").shape, (2,))
         assert_equal(np.percentile(x, (25, 50, 75),
-                                   interpolation="higher").shape, (3,))
+                                   method="higher").shape, (3,))
         assert_equal(np.percentile(x, (25, 50), axis=0,
-                                   interpolation="higher").shape, (2, 4, 5, 6))
+                                   method="higher").shape, (2, 4, 5, 6))
         assert_equal(np.percentile(x, (25, 50), axis=1,
-                                   interpolation="higher").shape, (2, 3, 5, 6))
+                                   method="higher").shape, (2, 3, 5, 6))
         assert_equal(np.percentile(x, (25, 50), axis=2,
-                                   interpolation="higher").shape, (2, 3, 4, 6))
+                                   method="higher").shape, (2, 3, 4, 6))
         assert_equal(np.percentile(x, (25, 50), axis=3,
-                                   interpolation="higher").shape, (2, 3, 4, 5))
+                                   method="higher").shape, (2, 3, 4, 5))
         assert_equal(np.percentile(x, (25, 50, 75), axis=1,
-                                   interpolation="higher").shape, (3, 3, 5, 6))
+                                   method="higher").shape, (3, 3, 5, 6))
 
     def test_scalar_q(self):
         # test for no empty dimensions for compatibility with old percentile
@@ -3057,33 +3057,33 @@ def test_scalar_q(self):
 
         # test for no empty dimensions for compatibility with old percentile
         x = np.arange(12).reshape(3, 4)
-        assert_equal(np.percentile(x, 50, interpolation='lower'), 5.)
+        assert_equal(np.percentile(x, 50, method='lower'), 5.)
         assert_(np.isscalar(np.percentile(x, 50)))
         r0 = np.array([4.,  5.,  6.,  7.])
-        c0 = np.percentile(x, 50, interpolation='lower', axis=0)
+        c0 = np.percentile(x, 50, method='lower', axis=0)
         assert_equal(c0, r0)
         assert_equal(c0.shape, r0.shape)
         r1 = np.array([1.,  5.,  9.])
-        c1 = np.percentile(x, 50, interpolation='lower', axis=1)
+        c1 = np.percentile(x, 50, method='lower', axis=1)
         assert_almost_equal(c1, r1)
         assert_equal(c1.shape, r1.shape)
 
         out = np.empty((), dtype=x.dtype)
-        c = np.percentile(x, 50, interpolation='lower', out=out)
+        c = np.percentile(x, 50, method='lower', out=out)
         assert_equal(c, 5)
         assert_equal(out, 5)
         out = np.empty(4, dtype=x.dtype)
-        c = np.percentile(x, 50, interpolation='lower', axis=0, out=out)
+        c = np.percentile(x, 50, method='lower', axis=0, out=out)
         assert_equal(c, r0)
         assert_equal(out, r0)
         out = np.empty(3, dtype=x.dtype)
-        c = np.percentile(x, 50, interpolation='lower', axis=1, out=out)
+        c = np.percentile(x, 50, method='lower', axis=1, out=out)
         assert_equal(c, r1)
         assert_equal(out, r1)
 
     def test_exception(self):
         assert_raises(ValueError, np.percentile, [1, 2], 56,
-                      interpolation='foobar')
+                      method='foobar')
         assert_raises(ValueError, np.percentile, [1], 101)
         assert_raises(ValueError, np.percentile, [1], -1)
         assert_raises(ValueError, np.percentile, [1], list(range(50)) + [101])
@@ -3124,12 +3124,12 @@ def test_percentile_out(self):
         # q.dim > 1, int
         r0 = np.array([[0,  1,  2, 3], [4, 5, 6, 7]])
         out = np.empty((2, 4), dtype=x.dtype)
-        c = np.percentile(x, (25, 50), interpolation='lower', axis=0, out=out)
+        c = np.percentile(x, (25, 50), method='lower', axis=0, out=out)
         assert_equal(c, r0)
         assert_equal(out, r0)
         r1 = np.array([[0,  4,  8], [1,  5,  9]])
         out = np.empty((2, 3), dtype=x.dtype)
-        c = np.percentile(x, (25, 50), interpolation='lower', axis=1, out=out)
+        c = np.percentile(x, (25, 50), method='lower', axis=1, out=out)
         assert_equal(c, r1)
         assert_equal(out, r1)
 
@@ -3146,10 +3146,10 @@ def test_percentile_empty_dim(self):
         assert_array_equal(np.percentile(d, 50, axis=-4).shape, (1, 2, 1))
 
         assert_array_equal(np.percentile(d, 50, axis=2,
-                                         interpolation='midpoint').shape,
+                                         method='midpoint').shape,
                            (11, 1, 1))
         assert_array_equal(np.percentile(d, 50, axis=-2,
-                                         interpolation='midpoint').shape,
+                                         method='midpoint').shape,
                            (11, 1, 1))
 
         assert_array_equal(np.array(np.percentile(d, [10, 50], axis=0)).shape,
@@ -3172,10 +3172,10 @@ def test_percentile_no_overwrite(self):
 
     def test_no_p_overwrite(self):
         p = np.linspace(0., 100., num=5)
-        np.percentile(np.arange(100.), p, interpolation="midpoint")
+        np.percentile(np.arange(100.), p, method="midpoint")
         assert_array_equal(p, np.linspace(0., 100., num=5))
         p = np.linspace(0., 100., num=5).tolist()
-        np.percentile(np.arange(100.), p, interpolation="midpoint")
+        np.percentile(np.arange(100.), p, method="midpoint")
         assert_array_equal(p, np.linspace(0., 100., num=5).tolist())
 
     def test_percentile_overwrite(self):
@@ -3253,14 +3253,14 @@ def test_out(self):
         o = np.zeros((4,))
         d = np.ones((3, 4))
         assert_equal(np.percentile(d, 0, 0, out=o), o)
-        assert_equal(np.percentile(d, 0, 0, interpolation='nearest', out=o), o)
+        assert_equal(np.percentile(d, 0, 0, method='nearest', out=o), o)
         o = np.zeros((3,))
         assert_equal(np.percentile(d, 1, 1, out=o), o)
-        assert_equal(np.percentile(d, 1, 1, interpolation='nearest', out=o), o)
+        assert_equal(np.percentile(d, 1, 1, method='nearest', out=o), o)
 
         o = np.zeros(())
         assert_equal(np.percentile(d, 2, out=o), o)
-        assert_equal(np.percentile(d, 2, interpolation='nearest', out=o), o)
+        assert_equal(np.percentile(d, 2, method='nearest', out=o), o)
 
     def test_out_nan(self):
         with warnings.catch_warnings(record=True):
@@ -3270,15 +3270,15 @@ def test_out_nan(self):
             d[2, 1] = np.nan
             assert_equal(np.percentile(d, 0, 0, out=o), o)
             assert_equal(
-                np.percentile(d, 0, 0, interpolation='nearest', out=o), o)
+                np.percentile(d, 0, 0, method='nearest', out=o), o)
             o = np.zeros((3,))
             assert_equal(np.percentile(d, 1, 1, out=o), o)
             assert_equal(
-                np.percentile(d, 1, 1, interpolation='nearest', out=o), o)
+                np.percentile(d, 1, 1, method='nearest', out=o), o)
             o = np.zeros(())
             assert_equal(np.percentile(d, 1, out=o), o)
             assert_equal(
-                np.percentile(d, 1, interpolation='nearest', out=o), o)
+                np.percentile(d, 1, method='nearest', out=o), o)
 
     def test_nan_behavior(self):
         a = np.arange(24, dtype=float)
@@ -3333,13 +3333,13 @@ def test_nan_behavior(self):
         b[:, 1] = np.nan
         b[:, 2] = np.nan
         assert_equal(np.percentile(a, [0.3, 0.6], (0, 2)), b)
-        # axis02 not zerod with nearest interpolation
+        # axis02 not zerod with method='nearest'
         b = np.percentile(np.arange(24, dtype=float).reshape(2, 3, 4),
-                          [0.3, 0.6], (0, 2), interpolation='nearest')
+                          [0.3, 0.6], (0, 2), method='nearest')
         b[:, 1] = np.nan
         b[:, 2] = np.nan
         assert_equal(np.percentile(
-            a, [0.3, 0.6], (0, 2), interpolation='nearest'), b)
+            a, [0.3, 0.6], (0, 2), method='nearest'), b)
 
     def test_nan_q(self):
         # GH18830
@@ -3412,26 +3412,36 @@ def test_no_p_overwrite(self):
         # this is worth retesting, because quantile does not make a copy
         p0 = np.array([0, 0.75, 0.25, 0.5, 1.0])
         p = p0.copy()
-        np.quantile(np.arange(100.), p, interpolation="midpoint")
+        np.quantile(np.arange(100.), p, method="midpoint")
         assert_array_equal(p, p0)
 
         p0 = p0.tolist()
         p = p.tolist()
-        np.quantile(np.arange(100.), p, interpolation="midpoint")
+        np.quantile(np.arange(100.), p, method="midpoint")
         assert_array_equal(p, p0)
 
     @pytest.mark.parametrize("dtype", np.typecodes["AllInteger"])
     def test_quantile_preserve_int_type(self, dtype):
         res = np.quantile(np.array([1, 2], dtype=dtype), [0.5],
-                          interpolation="nearest")
+                          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(
diff --git a/numpy/lib/tests/test_nanfunctions.py b/numpy/lib/tests/test_nanfunctions.py
index 126dba495a27..733a077ea1fc 100644
--- a/numpy/lib/tests/test_nanfunctions.py
+++ b/numpy/lib/tests/test_nanfunctions.py
@@ -1108,12 +1108,12 @@ def test_no_p_overwrite(self):
         # this is worth retesting, because quantile does not make a copy
         p0 = np.array([0, 0.75, 0.25, 0.5, 1.0])
         p = p0.copy()
-        np.nanquantile(np.arange(100.), p, interpolation="midpoint")
+        np.nanquantile(np.arange(100.), p, method="midpoint")
         assert_array_equal(p, p0)
 
         p0 = p0.tolist()
         p = p.tolist()
-        np.nanquantile(np.arange(100.), p, interpolation="midpoint")
+        np.nanquantile(np.arange(100.), p, method="midpoint")
         assert_array_equal(p, p0)
 
     @pytest.mark.parametrize("axis", [None, 0, 1])
diff --git a/numpy/typing/tests/data/reveal/lib_function_base.pyi b/numpy/typing/tests/data/reveal/lib_function_base.pyi
index 854b955b4a98..c559eb295893 100644
--- a/numpy/typing/tests/data/reveal/lib_function_base.pyi
+++ b/numpy/typing/tests/data/reveal/lib_function_base.pyi
@@ -144,7 +144,7 @@ reveal_type(np.percentile(AR_O, 50))  # E: Any
 reveal_type(np.percentile(AR_f8, [50]))  # E: ndarray[Any, dtype[floating[Any]]]
 reveal_type(np.percentile(AR_c16, [50]))  # E: ndarray[Any, dtype[complexfloating[Any, Any]]]
 reveal_type(np.percentile(AR_m, [50]))  # E: ndarray[Any, dtype[timedelta64]]
-reveal_type(np.percentile(AR_M, [50], interpolation="nearest"))  # E: ndarray[Any, dtype[datetime64]]
+reveal_type(np.percentile(AR_M, [50], method="nearest"))  # E: ndarray[Any, dtype[datetime64]]
 reveal_type(np.percentile(AR_O, [50]))  # E: ndarray[Any, dtype[object_]]
 reveal_type(np.percentile(AR_f8, [50], keepdims=True))  # E: Any
 reveal_type(np.percentile(AR_f8, [50], axis=[1]))  # E: Any
@@ -158,7 +158,7 @@ reveal_type(np.quantile(AR_O, 0.5))  # E: Any
 reveal_type(np.quantile(AR_f8, [0.5]))  # E: ndarray[Any, dtype[floating[Any]]]
 reveal_type(np.quantile(AR_c16, [0.5]))  # E: ndarray[Any, dtype[complexfloating[Any, Any]]]
 reveal_type(np.quantile(AR_m, [0.5]))  # E: ndarray[Any, dtype[timedelta64]]
-reveal_type(np.quantile(AR_M, [0.5], interpolation="nearest"))  # E: ndarray[Any, dtype[datetime64]]
+reveal_type(np.quantile(AR_M, [0.5], method="nearest"))  # E: ndarray[Any, dtype[datetime64]]
 reveal_type(np.quantile(AR_O, [0.5]))  # E: ndarray[Any, dtype[object_]]
 reveal_type(np.quantile(AR_f8, [0.5], keepdims=True))  # E: Any
 reveal_type(np.quantile(AR_f8, [0.5], axis=[1]))  # E: Any


Code review comments:
diff_hunk: @@ -1230,3 +1230,23 @@ def test_deprecated_module(self):
     def test_deprecated_attr(self):
         finfo = np.finfo(float)
         self.assert_deprecated(lambda: getattr(finfo, "machar"))
+
+
+class TestQuantileInterpolationDeprecation(_DeprecationTestCase):
+    # Deprecated 2021-11-08, NumPy 1.22
+    @pytest.mark.parametrize("func",
+        [np.percentile, np.quantile, np.nanpercentile, np.nanquantile])
+    def test_deprecated(self, func):
+        self.assert_deprecated(
+            lambda: func([0., 1.], 0., interpolation="linear"))
+        self.assert_deprecated(
+            lambda: func([0., 1.], 0., interpolation="nearest"))
+
+    @pytest.mark.parametrize("func",
+            [np.percentile, np.quantile, np.nanpercentile, np.nanquantile])
+    def test_both_passed(self, func):
+        with warnings.catch_warnings():
+            # catch the warning, but make sure it does not raise:
path: numpy/core/tests/test_deprecations.py
text: @user1:
"does not" <- "does"?

@author:
The warning is given before the error is thrown.  So we make sure the `DeprecationWarning` does _not_ raise, yes.  (In order to test the error that comes after the warning).

@user1:
Maybe "catch the DeprecationWarning  so that it does not raise".

diff_hunk: @@ -61,13 +61,13 @@
 # is made of a integer part (a.k.a 'i' or 'left') and a fractional part
 # (a.k.a 'g' or 'gamma')
 #
-# Each _QuantileInterpolation have two properties
+# Each _QuantileMethods has two properties
path: numpy/lib/function_base.py
text: @user1:
"Each method in"

@author:
```suggestion
# Each method in _QuantileMethods has two properties
```