Repository: scipy/scipy
Pull request: ENH: n-D and nan_policy support for scipy.stats.percentileofscore
ENH: Implement nd support for percentileofscore

The enhancement is inspired by https://stackoverflow.com/a/12414469,
which (as necessary before this) uses pure python loops, which is slow.

Note that BrenBarn's top-ranked answer on that question uses `rankdata`,
which has been added to the 'see also' section. The method of `rankdata`
uses `argsort` (as in an older version of `percentileofscore`) and so
it is faster. For ex., for assigning percentiles to each value in an
array of length n `rankdata` is `O(n log n)`,
while this `percentileofscore` is `O(n^2)`.
However, the algorithm of `rankdata` is not equivalent when
the list contains duplicates.

The enchantment also required an overhaul of the treatment of nan's.

Closes #2868

Problem statement:
nan and stats.percentileofscore
In function `stats.percentileofscore`, `nan` is (effectively) treated as +infinity.
E.g.

```
scipy.stats.percentileofscore ([3, 1, numpy.nan, numpy.nan, 2], 5)
```

returns `60.00` (as it the percent of values below 5).

However, other numerical functions (e.g. `scipy.mean`) return `nan` for entries containing `nan`s. 

Is it the desired behaviour of `stats.percentileofscore`?

Hints from issue discussion:
Functions that are not specifically designed for handling nans (*), don't check whether there are nans.
The behavior then is just a consequence of the implementation. In this case, the numpy `sort` function (which is used by `percentileofscore`, moves nans to the end. In (numpy.)`mean`, all numbers including nans are summed which results in a NaN as outcome.

(*) there are a few functions with nan in the name. However, for missing values you can use masked arrays and the scipy.stats.mstats functions (which unfortunately don't have currently the same signature and features as the "regular" stats functions.

So behavior with nans is "undefined" or arbitrary, not a desired behavior.

@josef-pkt `percentileofscore` (not be be confused with `scoreatpercentile`) does not use `sort`; it use comparisons of the form:

```
sum(a < score) / float(n) * 100
```

But anyway, thanks (I thought that there is some convention of dealing with fields with nans).

(I always confuse them and didn't check this time.)

sort is used if `kind == 'rank'`

inequality comparison with NaN are always False, so the outcome depends on the direction of the inequality and where the nan is

```
>>> 5 < np.array([3, 1, np.nan, np.nan, 2])
array([False, False, False, False, False], dtype=bool)
>>> 5 < np.nan
False
>>> 5 > np.nan
False
```

NaN checking is expensive, and I was arguing so far that we shouldn't do it by default. However, we might make NaN or missing value handling optional at some point in the future. (If someone is interested in working on this and we can make the masked array and stats versions consistent.)

On Fri, Sep 13, 2013 at 8:15 PM, Josef Perktold notifications@github.comwrote:

> Functions that are not specifically designed for handling nans (*), don't
> check whether there are nans.
> The behavior then is just a consequence of the implementation. In this
> case, the numpy sort function (which is used by percentileofscore, moves
> nans to the end. In (numpy.)mean, all numbers including nans are summed
> which results in a NaN as outcome.

Searchsorted also works with nan i.e.

In [1]: a = array([nan]_5 + [0]_5)

In [2]: a.sort()

In [3]: a
Out[3]: array([  0.,   0.,   0.,   0.,   0.,  nan,  nan,  nan,  nan,  nan])

In [4]: a.searchsorted(nan)
Out[4]: 5

Which might be useful for percentiles and such.

(*) there are a few functions with nan in the name. However, for missing

> values you can use masked arrays and the scipy.stats.mstats functions
> (which unfortunately don't have currently the same signature and features
> as the "regular" stats functions.
> 
> So behavior with nans is "undefined" or arbitrary, not a desired behavior.

Chuck

Commit message: Add a kwarg to switch to safer & heavier nan-check

This allows me to remove my near-duplicate function
inside percentileofscore().

Patch to review:
diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py
index 8df8ffa33e4d..d68cb471af72 100644
--- a/scipy/stats/stats.py
+++ b/scipy/stats/stats.py
@@ -216,16 +216,20 @@
            'page_trend_test', 'somersd']
 
 
-def _contains_nan(a, nan_policy='propagate'):
+def _contains_nan(a, nan_policy='propagate', summary_check=True):
     policies = ['propagate', 'raise', 'omit']
     if nan_policy not in policies:
         raise ValueError("nan_policy must be one of {%s}" %
                          ', '.join("'%s'" % s for s in policies))
     try:
-        # Calling np.sum to avoid creating a huge array into memory
-        # e.g. np.isnan(a).any()
-        with np.errstate(invalid='ignore'):
-            contains_nan = np.isnan(np.sum(a))
+        # The summary check avoids creating a (potentially huge) array.
+        # However, it will set contains_nan to True for (e.g.) [-inf, 0, +inf].
+        # If this is undesirabl, set summary_check to False instead.
+        if summary_check:
+            with np.errstate(invalid='ignore'):
+                contains_nan = np.isnan(np.sum(a))
+        else:
+            contains_nan = np.isnan(a).any()
     except TypeError:
         # This can happen when attempting to sum things which are not
         # numbers (e.g. as in the function `mode`). Try an alternative method:
@@ -1972,7 +1976,7 @@ def _compute_qth_percentile(sorted_, per, interpolation_method, axis):
     return np.add.reduce(sorted_[tuple(indexer)] * weights, axis=axis) / sumval
 
 
-def percentileofscore(a, score, kind='rank'):
+def percentileofscore(a, score, kind='rank', nan_policy='propagate'):
     """
     Compute the percentile rank of a score relative to a list of scores.
 
@@ -1983,9 +1987,9 @@ def percentileofscore(a, score, kind='rank'):
     Parameters
     ----------
     a : array_like
-        Array of scores to which `score` is compared.
-    score : int or float
-        Score that is compared to the elements in `a`.
+        Array to which `score` is compared.
+    score : array_like
+        Array to compute percentiles for.
     kind : {'rank', 'weak', 'strict', 'mean'}, optional
         Specifies the interpretation of the resulting score.
         The following options are available (default is 'rank'):
@@ -1999,6 +2003,13 @@ def percentileofscore(a, score, kind='rank'):
             strictly less than the given score are counted.
           * 'mean': The average of the "weak" and "strict" scores, often used
             in testing.  See https://en.wikipedia.org/wiki/Percentile_rank
+    nan_policy : {'propagate', 'raise', 'omit'}, optional
+        Defines how to handle when input contains nan.
+        The following options are available (default is 'propagate'):
+
+          * 'propagate': returns nan
+          * 'raise': throws an error
+          * 'omit': performs the calculations ignoring nan values
 
     Returns
     -------
@@ -2008,6 +2019,7 @@ def percentileofscore(a, score, kind='rank'):
     See Also
     --------
     numpy.percentile
+    scipy.stats.scoreatpercentile, scipy.stats.rankdata
 
     Examples
     --------
@@ -2038,28 +2050,90 @@ def percentileofscore(a, score, kind='rank'):
     >>> stats.percentileofscore([1, 2, 3, 3, 4], 3, kind='mean')
     60.0
 
+    Multidimensional score arrays are supported:
+
+    >>> stats.percentileofscore([1, 2, 3, 3, 4], [2, 3])
+    array([40., 70.])
+
+    The inputs can be infinite:
+
+    >>> stats.percentileofscore([-np.inf, 0, 1, np.inf], [1, 100])
+    array([75., 75.])
+
+    If `a` is empty, or containing only nan's (which are omitted),
+    the resulting percentiles are all nan:
+
+    >>> stats.percentileofscore([], 1)
+    nan
+    >>> stats.percentileofscore([np.nan, 1, 2], [1, 2], nan_policy='omit')
+    array([ 50., 100.])
+    >>> stats.percentileofscore([np.nan, np.nan], [1, 2], nan_policy='omit')
+    array([nan, nan])
+
+    The 'omit' policy does not make sense if there are nan's among the scores,
+    and so an error in that case.
     """
-    if np.isnan(score):
-        return np.nan
+
     a = np.asarray(a)
     n = len(a)
+    score = np.asarray(score)
+
+    # Nan treatment
+    cna, npa = _contains_nan(a, nan_policy, summary_check=False)
+    cns, nps = _contains_nan(score, nan_policy, summary_check=False)
+
+    if (cna or cns) and nan_policy == 'raise':
+        raise ValueError("The input contains nan values")
+
+    elif cns and nan_policy == "omit":
+        raise ValueError("The input scores contains nan values,"
+                         " which is incompatible with the 'omit' policy.")
+
+    # If a score is nan, then the output should be nan
+    if cns and nan_policy == "propagate":
+        score = ma.masked_where(np.isnan(score), score)
+
+    # Don't count nans
+    if cna and nan_policy == "omit":
+        a = ma.masked_where(np.isnan(a), a)
+        n = a.count()
+
+    # With `propagate`, an `a` with nans should always yield nan
+    if cna and nan_policy == "propagate":
+        n = 0
+
+    # Cannot compare to empty list ==> nan
     if n == 0:
-        return 100.0
-
-    if kind == 'rank':
-        left = np.count_nonzero(a < score)
-        right = np.count_nonzero(a <= score)
-        pct = (right + left + (1 if right > left else 0)) * 50.0/n
-        return pct
-    elif kind == 'strict':
-        return np.count_nonzero(a < score) / n * 100
-    elif kind == 'weak':
-        return np.count_nonzero(a <= score) / n * 100
-    elif kind == 'mean':
-        pct = (np.count_nonzero(a < score) + np.count_nonzero(a <= score)) / n * 50
-        return pct
+        perct = np.full_like(score, np.nan, dtype=float)
+
     else:
-        raise ValueError("kind can only be 'rank', 'strict', 'weak' or 'mean'")
+        # Prepare broadcasting
+        score = score[..., None]
+        count = lambda x: np.count_nonzero(x, -1)
+
+        # Compute
+        if kind == 'rank':
+            left  = count(a < score)
+            right = count(a <= score)
+            plus1 = left < right
+            perct = (left + right + plus1) * (50.0 / n)
+        elif kind == 'strict':
+            perct = count(a < score) * (100.0 / n)
+        elif kind == 'weak':
+            perct = count(a <= score) * (100.0 / n)
+        elif kind == 'mean':
+            left  = count(a < score)
+            right = count(a <= score)
+            perct = (left + right) * (50.0 / n)
+        else:
+            raise ValueError("kind can only be 'rank', 'strict', 'weak' or 'mean'")
+
+    # Re-insert nan's
+    perct = ma.filled(perct, np.nan)
+
+    if perct.ndim == 0:
+        return perct[()]
+    return perct
 
 
 HistogramResult = namedtuple('HistogramResult',
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index e73b2784ba7c..eb06593030e3 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -2829,6 +2829,32 @@ def test_percentileofscore():
 
     assert_raises(ValueError, pcos, [1, 2, 3, 3, 4], 3, kind='unrecognized')
 
+    # score *arrays*
+    a = 1 + np.arange(10)
+    assert_equal(pcos(a, [4, 5]), [40.0, 50.0])
+    assert_equal(pcos(a, [[4, 5], [6, 7]]), [[40.0, 50.0], [60.0, 70.0]])
+    for (kind, result) in [('mean',   [35.0, 45.0]),
+                           ('strict', [30.0, 40.0]),
+                           ('weak',   [40.0, 50.0])]:
+        assert_equal(pcos(a, [4, 5], kind=kind), result)
+
+    # infinite inputs
+    assert_equal(pcos([1,3],      [-np.inf, +np.inf]), [0.0, 100.0])
+    assert_equal(pcos([1,np.inf], [-np.inf, +np.inf]), [0.0, 100.0])
+    assert_equal(pcos([1,np.inf], [-np.inf, +np.inf], kind="strict"), [0.0, 50.0])
+    assert_equal(pcos([1,np.inf], [-np.inf, +np.inf], kind="mean"),   [0.0, 75.0])
+    assert_equal(pcos([1,np.inf], [-np.inf, +np.inf], kind="weak"),   [0.0, 100.0])
+
+    # nan's
+    assert_equal(pcos([], 1), np.nan)
+    assert_equal(pcos([np.nan, 1, 2],   [0, 1, 2], nan_policy='omit'),      [0, 50.0, 100.0])
+    assert_equal(pcos([np.nan, 2],      [0, 1, 2], nan_policy='propagate'), [np.nan, np.nan, np.nan])
+    assert_equal(pcos([np.nan, np.nan], [0, 1, 2], nan_policy='omit'),      [np.nan, np.nan, np.nan])
+    assert_equal(pcos([1, 2], [1, 2, np.nan], nan_policy='propagate'),      [50.0, 100.0, np.nan])
+    assert_raises(ValueError, pcos, [1, 2, 3, np.nan], [1, 2, 3], nan_policy='raise')
+    assert_raises(ValueError, pcos, [1, 2, 3], [1, 2, 3, np.nan], nan_policy='raise')
+    assert_raises(ValueError, pcos, [1, 2, 3, 3, 4], [np.nan, 2], nan_policy='omit')
+
 
 PowerDivCase = namedtuple('Case', ['f_obs', 'f_exp', 'ddof', 'axis',
                                    'chi2',     # Pearson's


Code review comments:
diff_hunk: @@ -2829,6 +2829,32 @@ def test_percentileofscore():
 
     assert_raises(ValueError, pcos, [1, 2, 3, 3, 4], 3, kind='unrecognized')
 
+    # score *arrays*
+    a = 1 + np.arange(10)
+    assert_equal(pcos(a, [4, 5]), [40.0, 50.0])
path: scipy/stats/tests/test_stats.py
text: @user2:
Would be better to use parametrization with `@user1.mark.parametrize` as most tests are similar.

@author:
Thanks for the input!

I followed the example of the other tests for percentileofscore.

I believe if I am to use `@user1.mark.parametrize` then I would have to move the tests out of `test_percentileofscore()`, which seems not so nice. Am I understanding you correctly?

@user2:
Yes this is not a problem to have multiple functions. You could even use a class to group all related percentile tests together (cleanest option IMO).

```python
class TestPercentileOfScore:
    def test_xxx(self):
        ...

    def test_warnings(self):
        ...
```

@author:
Ok, I've refactored the tests to use a class and `@user1.mark.parametrize`. Let me know what you think!

@user2:
It's better thanks! Also some cases were not covered which you fixed.

I compared both versions and if I am not mistaking it matches. I think I am just missing this one (I think the one called `test_large_numbers_missing` should have been it):

```python
    for kind in ('rank', 'mean', 'strict', 'weak'):
        assert_equal(
              pcos([10, 20, 30, 50, 60, 70, 80, 90, 100, 110],
                   40, kind=kind), 30.0)
```

A last thing, I know it is not good for presentation, but we should not align list elements. So make sure to have just a single space after a comma. Not less (the first lists), not more.

@author:
You're right, I had omitted that one (and wrongly named another one).

Also fixed the spacing, as requested.