ENH: Weighted quantile #9211

chunweiyuan · 2017-06-02T21:38:44Z

Creates quantile() function that allows weights input.

The following are plotted the same way as the graphs at the bottom of the np.percentile docs, but with an added weights argument:

weights=[1, 2, 3, 4]:

weights = [2, 1, 1, 2]:

numpy/lib/function_base.py

eric-wieser · 2017-06-03T00:12:40Z

numpy/lib/function_base.py

+    if weights is None:
+        Nx = ap.shape[-1]
+        indices = q * (Nx - 1)
+        indices = np.broadcast_to(indices, ap[..., 0].shape + q.shape)


Why not ap.shape[:-1]?

eric-wieser · 2017-06-03T00:14:02Z

numpy/lib/function_base.py

+        def _take1d(vec, inds):
+            return take(vec, inds)
+
+        vec_take = np.vectorize(_take1d, signature='(n),(m)->(m)')


Isn't this a superset of arraysort above? Just use this for both?

That might work, but I think having separate names (and signature) provides more clarity on what these two routines are trying to accomplish. I think It certainly helps understanding/debugging the code. Unless there's a real performance impact I prefer to keep it the way it is, though I'd defer to your judgment.

eric-wieser · 2017-06-03T00:16:28Z

numpy/lib/function_base.py

+        An array of weights associated with the values in `a`. Each value in
+        `a` contributes to the average according to its associated weight.
+        The weights array can either be 1-D (in which case its length must be
+        the size of `a` along the given axis) or of the same shape as `a`.


I think this should be relaxed to "Either 1D, or such that weights.ndim == a.ndim and that weights and a are broadcastable"

numpy/lib/function_base.py

juliantaylor · 2017-06-03T08:21:03Z

numpy/lib/function_base.py

+        indices = q * (Nx - 1)
+        indices = np.broadcast_to(indices, ap.shape[:-1] + q.shape)
+        ind_sorted = np.argsort(ap, axis=-1)  # sort values long axis
+        ap = arraysort(ap, ind_sorted)


without weights it should not be slower than the old percentile, so partition should be used and the vectorization of the index sort avoided.

Just fixed it in the latest push.

juliantaylor · 2017-06-03T08:27:27Z

this does seem almost out of scope for numpy and more something for a statistics package, has this been previously discussed?

bashtage · 2017-06-03T11:42:08Z

The explanation of weights can be substantially simplified if you make the connection between the normalized weights and repeating the data in the sample in the integer case.

For example, if you have [1,2,3] with weight [1,2,3] this will be have identically to the case where there are no weighs and the data are [1,2,2,3,3,3]. In Stata, for example, this type if weighting is known as fweighting (for frequency). There are other interpretations of weights that may result in a different estimator.

chunweiyuan · 2017-06-03T21:51:27Z

I think this PR belongs in numpy:
1.) it's nothing more than an extension of percentile, and
2.) I've gone through great pain to make sure the underlying interpolation scheme, when weights is not None, is consistent with the old percentile. In the new tests there are cases where all weights are identical, and quantile renders the same result as the old percentile. There's another case where weights=[0, 1, 2], and quantile gives the same result as using the percentile on a simple repetition of the original data.

I've tried to address the comments above in the latest push. The explanation of how I compute the weighted quantile is still a bit cryptic, but I can't find a better way at the moment.

chunweiyuan · 2017-07-20T20:22:29Z

Anything I could do to help expedite this? It seems that the merge conflict in #9213 is rather small, and I could help bridge it with this PR as well.

eric-wieser · 2017-07-20T20:25:16Z

8000

@chunweiyuan: Thanks for pointing out the merge conflict, I can go fix that.

There was a lack of decision in the comments of #9213 whether quantile should be in the API, if we already have percentile. Perhaps asking the mailing list for opinions would resolve that.

chunweiyuan · 2017-08-03T17:10:41Z

Doesn't seem like there's much response to my email related to the API, so I think the way forward is to keep both quantile and percentile. There was, however, one suggestion related to using an O(n) partition algorithm instead of O(n log n) (argsort) in the code. I'm not sure how that would be addressed since I need ind_sorted in my code to sort the weights array.

madphysicist · 2017-08-04T18:15:37Z

numpy/lib/function_base.py

@@ -4202,61 +4239,166 @@ def percentile(a, q, axis=None, out=None,
    Returns
    -------
    percentile : scalar or ndarray


This should be quantile, not percentile

Thanks. Fixed. This baby is now just waiting for #9213........

madphysicist

This probably needs an entry somewhere in the release notes.

madphysicist · 2017-08-09T18:45:10Z

numpy/lib/tests/test_function_base.py

@@ -3185,6 +3185,85 @@ def test_nan_behavior(self):
                a, [0.3, 0.6], (0, 2), interpolation='nearest'), b)


+class TestQuantile(object):
+
+    def test_reqression(self):


This looks like a q instaead of a g. If I am wrong, ignore.

You're right.

madphysicist · 2017-08-09T18:45:52Z

numpy/lib/tests/test_function_base.py

+        actual = np.quantile(ar, q=0.5, weights=np.ones(8).reshape(2, 1, 4))
+        assert_almost_equal(actual, expected)
+
+    def test_various_weights(self):


How about integer weights?

I think line 3233 there tests integer weights.

madphysicist · 2017-08-09T18:50:39Z

numpy/lib/function_base.py

    """
    Compute the qth percentile of the data along the specified axis.

    Returns the qth percentile(s) of the array elements.

+    Uses quantile() underneath, with
+    1.) weights=None, and


Remove the parens.

madphysicist · 2017-08-09T18:55:29Z

numpy/lib/function_base.py

    """
    Compute the qth percentile of the data along the specified axis.

    Returns the qth percentile(s) of the array elements.

+    Uses quantile() underneath, with


This should be removed or modified. There does not need to be any implementation information in the docs. Instead, you could say something like "Calling percentile is identical to calling quantile with weights set to None and the percentiles divided by 100."

madphysicist · 2017-08-09T18:58:12Z

numpy/lib/function_base.py

+    --------
+    quantile
+    """
+    return quantile(a, q=np.asarray(q)/100., axis=axis, weights=None,


You could save on a call to asarray by doing q=true_divide(q, 100)

madphysicist · 2017-08-10T21:24:45Z

@chunweiyuan For future reference, it is probably better to rebase your branch onto master rather than merging master into your branch. The history will look much cleaner and you won't have to deal with a stray merge commit later down the line.

chunweiyuan · 2017-08-10T21:47:10Z

Sounds good. Old habits die hard. Now this branch is just waiting on #9213, I believe.

eric-wieser · 2017-08-11T18:53:33Z

Looks like you did something screwy with a merge/rebase there. I've seen that behaviour before, but don't remember the cause.

eric-wieser · 2017-08-11T18:55:10Z

Ok, looks like actually github was being screwy - changing the base twice has removed the commits from #9475 from this PR, which was what was weird in the first place.

…n lib.function_base and nanquantile/nanpercentile in lib.nanfunctions (numpy#9211). Move weights normalization to function_base._validate_and_ureduce_weights(), and speed improvement to conditional statement therein (numpy#9211). Additional docs to how weights are applied for weighted quantile calculation in lib.function_base (numpy#9211). Consolidate all np.vectorize() algo into one in lib.function_base._quantile() (numpy#9211). Better mapping from weight space previous/next indexes to real space in lib.function_base._quantile(), with additional comments. (numpy#9211)