diff options
-rw-r--r-- | doc/release/1.9.0-notes.rst | 2 | ||||
-rw-r--r-- | numpy/lib/function_base.py | 161 | ||||
-rw-r--r-- | numpy/lib/tests/test_function_base.py | 113 |
3 files changed, 192 insertions, 84 deletions
diff --git a/doc/release/1.9.0-notes.rst b/doc/release/1.9.0-notes.rst index b78211d24..a958fe705 100644 --- a/doc/release/1.9.0-notes.rst +++ b/doc/release/1.9.0-notes.rst @@ -19,6 +19,8 @@ Future Changes Compatibility notes =================== +numpy.percentile returns an array instead of a list. + New Features ============ diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index ada54135e..472d7eecc 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -16,7 +16,8 @@ import sys import numpy.core.numeric as _nx from numpy.core import linspace from numpy.core.numeric import ones, zeros, arange, concatenate, array, \ - asarray, asanyarray, empty, empty_like, ndarray, around + asarray, asanyarray, empty, empty_like, ndarray, around, floor, \ + ceil, take from numpy.core.numeric import ScalarType, dot, where, newaxis, intp, \ integer, isscalar from numpy.core.umath import pi, multiply, add, arctan2, \ @@ -2765,7 +2766,9 @@ def median(a, axis=None, out=None, overwrite_input=False): # and check, use out array. return mean(part[indexer], axis=axis, out=out) -def percentile(a, q, axis=None, out=None, overwrite_input=False): + +def percentile(a, q, interpolation='linear', axis=None, out=None, + overwrite_input=False): """ Compute the qth percentile of the data along the specified axis. @@ -2777,29 +2780,40 @@ def percentile(a, q, axis=None, out=None, overwrite_input=False): Input array or object that can be converted to an array. q : float in range of [0,100] (or sequence of floats) Percentile to compute which must be between 0 and 100 inclusive. + interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'} + This optional parameter specifies the interpolation method to use, + when the desired quantile lies between two data points `i` and `j`: + * linear: `i + (j - i) * fraction`, where `fraction` is the + fractional part of the index surrounded by `i` and `j`. + * lower: `i`. + * higher: `j`. + * nearest: `i` or `j` whichever is nearest. + * midpoint: (`i` + `j`) / 2. axis : int, optional Axis along which the percentiles are computed. The default (None) - is to compute the median along a flattened version of the array. + is to compute the percentiles along a flattened version of the array. out : ndarray, optional Alternative output array in which to place the result. It must have the same shape and buffer length as the expected output, but the type (of the output) will be cast if necessary. overwrite_input : bool, optional - If True, then allow use of memory of input array `a` for - calculations. The input array will be modified by the call to - median. This will save memory when you do not need to preserve - the contents of the input array. Treat the input as undefined, - but it will probably be fully or partially sorted. - Default is False. Note that, if `overwrite_input` is True and the - input is not already an array, an error will be raised. + If True, then allow use of memory of input array `a` for + calculations. The input array will be modified by the call to + percentile. This will save memory when you do not need to preserve + the contents of the input array. In this case you should not make + any assumptions about the content of the passed in array `a` after + this function completes -- treat it as undefined. Default is False. + Note that, if the `a` input is not already an array this parameter + will have no effect, `a` will be converted to an array internally + regardless of the value of this parameter. Returns ------- - pcntile : ndarray + percentile : ndarray A new array holding the result (unless `out` is specified, in - which case that array is returned instead). If the input contains + which case that array is returned instead). If the input contains integers, or floats of smaller precision than 64, then the output - data-type is float64. Otherwise, the output data-type is the same + data-type is float64. Otherwise, the output data-type is the same as that of the input. See Also @@ -2809,10 +2823,11 @@ def percentile(a, q, axis=None, out=None, overwrite_input=False): Notes ----- Given a vector V of length N, the qth percentile of V is the qth ranked - value in a sorted copy of V. A weighted average of the two nearest - neighbors is used if the normalized ranking does not match q exactly. - 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``. + value in a sorted copy of V. The values and distances of the two nearest + neighbors as well as the `interpolation` parameter will determine the + percentile if the normalized ranking does not match 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``. Examples -------- @@ -2821,84 +2836,94 @@ def percentile(a, q, axis=None, out=None, overwrite_input=False): array([[10, 7, 4], [ 3, 2, 1]]) >>> np.percentile(a, 50) - 3.5 + array([ 3.5]) >>> np.percentile(a, 50, axis=0) - array([ 6.5, 4.5, 2.5]) + array([[ 6.5, 4.5, 2.5]]) >>> np.percentile(a, 50, axis=1) - array([ 7., 2.]) + array([[ 7.], + [ 2.]]) >>> m = np.percentile(a, 50, axis=0) >>> out = np.zeros_like(m) >>> np.percentile(a, 50, axis=0, out=m) - array([ 6.5, 4.5, 2.5]) + array([[ 6.5, 4.5, 2.5]]) >>> m - array([ 6.5, 4.5, 2.5]) + array([[ 6.5, 4.5, 2.5]]) >>> b = a.copy() >>> np.percentile(b, 50, axis=1, overwrite_input=True) - array([ 7., 2.]) + array([[ 7.], + [ 2.]]) >>> assert not np.all(a==b) >>> b = a.copy() >>> np.percentile(b, 50, axis=None, overwrite_input=True) - 3.5 + array([ 3.5]) """ - a = np.asarray(a) - - if q == 0: - return a.min(axis=axis, out=out) - elif q == 100: - return a.max(axis=axis, out=out) + a = asarray(a) + q = atleast_1d(q) + q = q / 100.0 + if (q < 0).any() or (q > 1).any(): + raise ValueError("Percentiles must be in the range [0,100]") + # prepare a for partioning if overwrite_input: if axis is None: - sorted = a.ravel() - sorted.sort() + ap = a.ravel() else: - a.sort(axis=axis) - sorted = a + ap = a else: - sorted = sort(a, axis=axis) + if axis is None: + ap = a.flatten() + else: + ap = a.copy() + if axis is None: axis = 0 - return _compute_qth_percentile(sorted, q, axis, out) + Nx = ap.shape[axis] + indices = q * (Nx - 1) + + # round fractional indices according to interpolation method + if interpolation == 'lower': + indices = floor(indices).astype(intp) + elif interpolation == 'higher': + indices = ceil(indices).astype(intp) + elif interpolation == 'midpoint': + indices = floor(indices) + 0.5 + elif interpolation == 'nearest': + indices = around(indices).astype(intp) + elif interpolation == 'linear': + pass # keep index as fraction and interpolate + else: + raise ValueError("interpolation can only be 'linear', 'lower' " + "'higher', 'midpoint', or 'nearest'") + + if indices.dtype == intp: # take the points along axis + ap.partition(indices, axis=axis) + return take(ap, indices, axis=axis, out=out) + else: # weight the points above and below the indices + indices_below = floor(indices).astype(intp) + indices_above = indices_below + 1 + indices_above[indices_above > Nx - 1] = Nx - 1 -# handle sequence of q's without calling sort multiple times -def _compute_qth_percentile(sorted, q, axis, out): - if not isscalar(q): - p = [_compute_qth_percentile(sorted, qi, axis, None) - for qi in q] + weights_above = indices - indices_below + weights_below = 1.0 - weights_above - if out is not None: - out.flat = p + weights_shape = [1, ] * ap.ndim + weights_shape[axis] = len(indices) + weights_below.shape = weights_shape + weights_above.shape = weights_shape - return p + ap.partition(concatenate((indices_below, indices_above)), axis=axis) + x1 = take(ap, indices_below, axis=axis) * weights_below + x2 = take(ap, indices_above, axis=axis) * weights_above + + if out is not None: + return add(x1, x2, out=out) + else: + return add(x1, x2) - q = q / 100.0 - if (q < 0) or (q > 1): - raise ValueError("percentile must be either in the range [0,100]") - - indexer = [slice(None)] * sorted.ndim - Nx = sorted.shape[axis] - index = q*(Nx-1) - i = int(index) - if i == index: - indexer[axis] = slice(i, i+1) - weights = array(1) - sumval = 1.0 - else: - indexer[axis] = slice(i, i+2) - j = i + 1 - weights = array([(j - index), (index - i)], float) - wshape = [1]*sorted.ndim - wshape[axis] = 2 - weights.shape = wshape - sumval = weights.sum() - - # Use add.reduce in both cases to coerce data type as well as - # check and use out array. - return add.reduce(sorted[indexer]*weights, axis=axis, out=out)/sumval def trapz(y, x=None, dx=1.0, axis=-1): """ diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index f52eb5fbe..dd0b6e0ee 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1440,27 +1440,108 @@ def compare_results(res, desired): assert_array_equal(res[i], desired[i]) -def test_percentile_list(): - assert_equal(np.percentile([1, 2, 3], 0), 1) +class TestScoreatpercentile(TestCase): + def test_basic(self): + x = np.arange(8) * 0.5 + assert_equal(np.percentile(x, 0), 0.) + assert_equal(np.percentile(x, 100), 3.5) + assert_equal(np.percentile(x, 50), 1.75) + + def test_2D(self): + x = np.array([[1, 1, 1], + [1, 1, 1], + [4, 4, 3], + [1, 1, 1], + [1, 1, 1]]) + assert_array_equal(np.percentile(x, 50, axis=0), [[1, 1, 1]]) + + def test_linear(self): + + # Test defaults + assert_equal(np.percentile(range(10), 50), 4.5) + + # explicitly specify interpolation_method 'fraction' (the default) + assert_equal(np.percentile(range(10), 50, + interpolation='linear'), 4.5) + + def test_lower_higher(self): + + # interpolation_method 'lower'/'higher' + assert_equal(np.percentile(range(10), 50, + interpolation='lower'), 4) + assert_equal(np.percentile(range(10), 50, + interpolation='higher'), 5) + + def test_midpoint(self): + assert_equal(np.percentile(range(10), 51, + interpolation='midpoint'), 4.5) + + def test_nearest(self): + assert_equal(np.percentile(range(10), 51, + interpolation='nearest'), 5) + assert_equal(np.percentile(range(10), 49, + interpolation='nearest'), 4) + + def test_sequence(self): + x = np.arange(8) * 0.5 + assert_equal(np.percentile(x, [0, 100, 50]), [0, 3.5, 1.75]) + + def test_axis(self): + x = np.arange(12).reshape(3, 4) + + assert_equal(np.percentile(x, (25, 50, 100)), [2.75, 5.5, 11.0]) + + r0 = [[2, 3, 4, 5], [4, 5, 6, 7], [8, 9, 10, 11]] + assert_equal(np.percentile(x, (25, 50, 100), axis=0), r0) + + r1 = [[0.75, 1.5, 3], [4.75, 5.5, 7], [8.75, 9.5, 11]] + assert_equal(np.percentile(x, (25, 50, 100), axis=1), r1) + + def test_exception(self): + assert_raises(ValueError, np.percentile, [1, 2], 56, + interpolation='foobar') + assert_raises(ValueError, np.percentile, [1], 101) + assert_raises(ValueError, np.percentile, [1], -1) + + def test_percentile_list(self): + assert_equal(np.percentile([1, 2, 3], 0), 1) + + def test_percentile_out(self): + x = np.array([1, 2, 3]) + y = np.zeros((3,)) + p = (1, 2, 3) + np.percentile(x, p, out=y) + assert_equal(y, np.percentile(x, p)) + + x = np.array([[1, 2, 3], + [4, 5, 6]]) + + y = np.zeros((3, 3)) + np.percentile(x, p, axis=0, out=y) + assert_equal(y, np.percentile(x, p, axis=0)) + + y = np.zeros((2, 3)) + np.percentile(x, p, axis=1, out=y) + assert_equal(y, np.percentile(x, p, axis=1)) + + def test_percentile_no_overwrite(self): + a = np.array([2, 3, 4, 1]) + np.percentile(a, [50], overwrite_input=False) + assert_equal(a, np.array([2, 3, 4, 1])) -def test_percentile_out(): - x = np.array([1, 2, 3]) - y = np.zeros((3,)) - p = (1, 2, 3) - np.percentile(x, p, out=y) - assert_equal(y, np.percentile(x, p)) + a = np.array([2, 3, 4, 1]) + np.percentile(a, [50]) + assert_equal(a, np.array([2, 3, 4, 1])) - x = np.array([[1, 2, 3], - [4, 5, 6]]) + def test_percentile_overwrite(self): + a = np.array([2, 3, 4, 1]) + b = np.percentile(a, [50], overwrite_input=True) + assert_equal(b, np.array([2.5])) - y = np.zeros((3, 3)) - np.percentile(x, p, axis=0, out=y) - assert_equal(y, np.percentile(x, p, axis=0)) + b = np.percentile([2, 3, 4, 1], [50], overwrite_input=True) + assert_equal(b, np.array([2.5])) - y = np.zeros((3, 2)) - np.percentile(x, p, axis=1, out=y) - assert_equal(y, np.percentile(x, p, axis=1)) class TestMedian(TestCase): |