Repository: scipy/scipy
Pull request: Add location and sign to KS test result
#### Reference issue
Aims to close #9389

#### What does this implement/fix?
Moves the result of Kolmogorov-Smirnov tests (1 or 2 samples) to a "tuple bunch" with two additional attributes:
- `statistic_location` indicating at what value of the stochastic variable (approximately) the maximal distance in CDF is reached
- `statistic_sign`, which is only nonempty for `alternative="two-sided"`, indicating whether the first distribution (for 2samp)/the empirical distribution (for 1samp) was above or below the second/theoretical distribution at the `statistic_location`. This is a cheap way to determine not only whether the distros are different, but also if they are, in what direction.

#### Additional information
Matt and I quickly discussed part 1, the `statistic_location`. I later realized having part 2 i.e. the sign is equally important for me and perhaps for other folks, and is a cheap addition in terms of performance.

I haven't adapted the part of the docstring that represents doctests yet, happy to do so but perhaps it'd be useful if one of the maintainers took a quick look the the draft PR first so I know if there are any major issues.

Thank you for your time and patience


Problem statement:
Kolmogorov Smirnov 2 samples returning max distance location [Enhancement]
<!-- 
Thank you for taking the time to report a SciPy issue.
Please describe the issue in detail, and for bug reports
fill in the fields below. You can delete the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
-->

2 sample Kolmogorov Smirnov test (KS test) is about finding the maximal vertical distance between two empirical cumulative distributions; a large max distance is unlikely under the null assumption that the two empirical distributions are sampled from the same theoretical curve. This is implemented in `scipy.stats.ks_twosamp`.

The function calculates but does not return the location of the cumulative at which the max distance happens, but that would be extremely useful for a number of machine learning algorithms, e.g. to get an initial guess for a classifier between the two kinds of data.

I propose to add an optional argument to the function to return the location of the max distance between the distributions in terms of indices of the data and also value of the quantitative variable.

### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be marked as a code block automatically 
-->
Here a rough initial implementation for example purposes. The final few lines are the only ones that differ from the current implementation.
```python
def max_distance_cumulative(data1, data2):
    """
    Computes the maximal vertical distance between cumulative distributions
    (this is the statistic for KS tests). Code mostly copied from
    scipy.stats.ks_twosamp

    Parameters
    ----------
    data1 : array_like
        First data set
    data2 : array_like
        Second data set
    Returns
    -------
    d : float
        Max distance, i.e. value of the Kolmogorov Smirnov test. Sign is + if
        the cumulative of data1 < the one of data2 at that location, else -.
    x : float
        Value of x where maximal distance d is reached.
    """
    from numpy import ma

    (data1, data2) = (ma.asarray(data1), ma.asarray(data2))
    (n1, n2) = (data1.count(), data2.count())
    mix = ma.concatenate((data1.compressed(), data2.compressed()))
    mixsort = mix.argsort(kind='mergesort')
    csum = np.where(mixsort < n1, 1./n1, -1./n2).cumsum()

    # Check for ties
    if len(np.unique(mix)) < (n1+n2):
        ind = np.r_[np.diff(mix[mixsort]).nonzero()[0], -1]
        csum = csum[ind]
        mixsort = mixsort[ind]

    csumabs = ma.abs(csum)
    i = csumabs.argmax()

    d = csum[i]
    # mixsort[i] contains the index of mix with the max distance
    x = mix[mixsort[i]]

    return (d, x)

```

I'm happy to make a PR if there is interest for this.

Hints from issue discussion:
I thought this suggestion is very useful to know the index at maxima.
If this idea put it in `scipy.stats.mstats.ks_twosamp()`, the return variable will be increased by 3(ks value, p-value, index at maxima).

If you do not mind, I would like to take this issue.
While I see the point that this information can be helpful in certain cases, I would not add it:

- it does not belong to the KS-test, i.e. most users don't need this Information when running the test (i meant in most cases when users just want to test whether the distributions are the same)
- the location is often not unique
- it is neither hard nor expensive to get all the locations (in the end, it is just argmax of cdf1 - cdf2) outside of the KS-test

Thank you Chris for your comment that revived the issue. While I respect your point of view, I don't think those points invalidate the issue:

 - it does belong to the test just like the statistics and p value, this location is used to calculate the test statistics directly. Moreover, although many users won't need it, there are certainly specific applications such as model building that are quite major

- non uniqueness is ok, we can report one just like argmax does

 - it is expensive in terms of expertise and boilerplate, since one has to take the scipy function, figure out what's going on (docs are not great there I'm afraid), and back engineer it. I've got a master in physics and a PhD in bioinformatics and it took me decent effort to retrace the steps, I don't think this is a trivial issue.

In my opinion this issue is a typical case where we could just expand the usefulness of the codebase with minimal effort. These situations are all too rare and we should just grab the opportunity when we can.
There is also the problem of backwards-compatability if a new element is added.

This is the code in `ks_2samp` to get the distance between the cdfs:
```
data1 = np.sort(data1)
data2 = np.sort(data2)
n1 = data1.shape[0]
n2 = data2.shape[0]
data_all = np.concatenate([data1, data2])
cdf1 = np.searchsorted(data1, data_all, side='right') / n1
cdf2 = np.searchsorted(data2, data_all, side='right') / n2
d = np.max(np.absolute(cdf1 - cdf2))
```

Let's see if there are further comments from other developers
gh-9753 changes lots of the 2-sample KS code, but only the last line in @chrisb83's snippet above is different.  Even so that PR could compute the argmin, argmax (for D+, D-) without much extra effort.

What would this mean in the 1-sample case? [Context: I've been thinking about extending `kstest` to cover `ks_2samp`, so that one could just call `kstest(data0, data1,...)` where`data0`, `data1` can be one of `(rvs0, cdf), (cdf0, cdf), (rvs0, rvs1)`. It already supports the first two.]



@iosonofabio Thanks for your code. It is exactly what I was looking for. However, I have a question about the code. I am feeding two CDF as numpy to the code, but I keep getting an error of " 'numpy.ndarray' object has no attribute 'count'". I have tried data1['Values'].count() as well. However, I still give the same error. I wonder if you or anyone else who had the same issue can guide me on how to get rid of this error.


Thank you so much.
This issue has been quiet for a while but I was exactly looking for this kind of functionality in `ks_2samp`. I am a modeller and knowing the location of the maximum is important in the context of my work (as @iosonofabio also suggested above). I will use the code above and would certainly appreciate this PR. 
I'm still happy to make a PR for this if there is hope. And since we are on the topic, it would seem reasonable to add and option to `ks_2samp` that returns the KS statistic with a sign instead of the absolute value. Would be SO USEFUL in many occasions.
given that there is interest in this functionality: instead of changing ks_2samp, what about adding a new function, e.g. ks_dist that returns the distance and the location?
Sure, I could do that, in fact I have the function already...

On Fri, Jul 10, 2020, at 06:33, Christoph Baumgarten wrote:
> 
> 
> given that there is interest in this functionality: instead of changing 
> ks_2samp, what about adding a new function, e.g. ks_dist that returns 
> the distance and the location?
> 
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub 
> <https://github.com/scipy/scipy/issues/9389#issuecomment-656337235>, or 
> unsubscribe 
> <https://github.com/notifications/unsubscribe-auth/AAJFEADVHYHWA6GSOJADXN3R2YSQ5ANCNFSM4F4T4CVA>.
>

Commit message: Add location and sign to KS test result

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index 6446a3140a9b..159496b875cd 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -7549,7 +7549,8 @@ def chisquare(f_obs, f_exp=None, ddof=0, axis=0):
                             lambda_="pearson")
 
 
-KstestResult = namedtuple('KstestResult', ('statistic', 'pvalue'))
+KstestResult = _make_tuple_bunch(['statistic', 'pvalue'],
+                                 ['statistic_location', 'statistic_sign'])
 
 
 def _compute_dplus(cdfvals):
@@ -7565,7 +7566,10 @@ def _compute_dplus(cdfvals):
       Maximum distance of the CDF values below Uniform(0, 1)
 """
     n = len(cdfvals)
-    return (np.arange(1.0, n + 1) / n - cdfvals).max()
+    dplus = (np.arange(1.0, n + 1) / n - cdfvals)
+    amax = dplus.argmax()
+    loc_max = cdfvals[amax]
+    return (dplus[amax], loc_max)
 
 
 def _compute_dminus(cdfvals):
@@ -7582,7 +7586,10 @@ def _compute_dminus(cdfvals):
 
     """
     n = len(cdfvals)
-    return (cdfvals - np.arange(0.0, n)/n).max()
+    dminus = (cdfvals - np.arange(0.0, n)/n)
+    amax = dminus.argmax()
+    loc_max = cdfvals[amax]
+    return (dminus[amax], loc_max)
 
 
 @_rename_parameter("mode", "method")
@@ -7617,11 +7624,23 @@ def ks_1samp(x, cdf, args=(), alternative='two-sided', method='auto'):
 
     Returns
     -------
-    statistic : float
-        KS test statistic, either D, D+ or D- (depending on the value
-        of 'alternative')
-    pvalue : float
-        One-tailed or two-tailed p-value.
+    res: KstestResult
+        An object containing attributes:
+
+        statistic : float
+            KS test statistic, either D, D+ or D- (depending on the value
+            of 'alternative')
+        pvalue : float
+            One-tailed or two-tailed p-value.
+        statistic_location: float
+            Value of the stochastic variable closest to the KS statistic. By
+            definition, the distance between the theoretical and sampled
+            distributions at this value is largest.
+        statistic_sign: str
+            For alternative == 'two-sided', this attribute is set to '-' if the
+            sampled distribution is above the theoretical one at the statistic
+            location, otherwise '+'. For other alternative values, it is an
+            empty string.
 
     See Also
     --------
@@ -7705,18 +7724,32 @@ def ks_1samp(x, cdf, args=(), alternative='two-sided', method='auto'):
     x = np.sort(x)
     cdfvals = cdf(x, *args)
 
+    # TODO: add logic to find the argmax
+
     if alternative == 'greater':
-        Dplus = _compute_dplus(cdfvals)
-        return KstestResult(Dplus, distributions.ksone.sf(Dplus, N))
+        Dplus, d_location = _compute_dplus(cdfvals)
+        return KstestResult(Dplus, distributions.ksone.sf(Dplus, N),
+                            statistic_location=d_location,
+                            statistic_sign='')
 
     if alternative == 'less':
-        Dminus = _compute_dminus(cdfvals)
-        return KstestResult(Dminus, distributions.ksone.sf(Dminus, N))
+        Dminus, d_location = _compute_dminus(cdfvals)
+        return KstestResult(Dminus, distributions.ksone.sf(Dminus, N),
+                            statistic_location=d_location,
+                            statistic_sign='')
 
     # alternative == 'two-sided':
-    Dplus = _compute_dplus(cdfvals)
-    Dminus = _compute_dminus(cdfvals)
-    D = np.max([Dplus, Dminus])
+    Dplus, dplus_location = _compute_dplus(cdfvals)
+    Dminus, dminus_location = _compute_dminus(cdfvals)
+    if Dplus > Dminus:
+        D = Dplus
+        d_location = dplus_location
+        d_sign = '+'
+    else:
+        D = Dminus
+        d_location = dminus_location
+        d_sign = '-'
+
     if mode == 'auto':  # Always select exact
         mode = 'exact'
     if mode == 'exact':
@@ -7727,7 +7760,9 @@ def ks_1samp(x, cdf, args=(), alternative='two-sided', method='auto'):
         # mode == 'approx'
         prob = 2 * distributions.ksone.sf(D, N)
     prob = np.clip(prob, 0, 1)
-    return KstestResult(D, prob)
+    return KstestResult(D, prob,
+                        statistic_location=d_location,
+                        statistic_sign=d_sign)
 
 
 Ks_2sampResult = KstestResult
@@ -7920,10 +7955,20 @@ def ks_2samp(data1, data2, alternative='two-sided', method='auto'):
 
     Returns
     -------
-    statistic : float
-        KS statistic.
-    pvalue : float
-        One-tailed or two-tailed p-value.
+    res: KstestResult
+        An object containing attributes:
+
+        statistic : float
+            KS test statistic.
+        pvalue : float
+            One-tailed or two-tailed p-value.
+        statistic_location: float
+            Value of the stochastic variable closest to the KS statistic. By
+            definition, the distance between the two sampled distributions at
+            this value is largest.
+        statistic_sign: str
+            '-' if the first distribution is above the second at the statistic
+            location, otherwise '+'.
 
     See Also
     --------
@@ -8054,11 +8099,35 @@ def ks_2samp(data1, data2, alternative='two-sided', method='auto'):
     cdf1 = np.searchsorted(data1, data_all, side='right') / n1
     cdf2 = np.searchsorted(data2, data_all, side='right') / n2
     cddiffs = cdf1 - cdf2
+
+    # Identify the location of the statistic
+    argminS = np.argmin(cddiffs)
+    argmaxS = np.argmax(cddiffs)
+    loc_minS = data_all[argminS]
+    loc_maxS = data_all[argmaxS]
+
     # Ensure sign of minS is not negative.
-    minS = np.clip(-np.min(cddiffs), 0, 1)
-    maxS = np.max(cddiffs)
-    alt2Dvalue = {'less': minS, 'greater': maxS, 'two-sided': max(minS, maxS)}
-    d = alt2Dvalue[alternative]
+    minS = np.clip(-cddiffs[argminS], 0, 1)
+    maxS = cddiffs[argmaxS]
+
+    if alternative == 'less':
+        d = minS
+        d_location = loc_minS
+        d_sign = '-'
+    elif alternative == 'greater':
+        d = maxS
+        d_location = loc_maxS
+        d_sign = '+'
+    else:
+        if minS > maxS:
+            d = minS
+            d_location = loc_minS
+            d_sign = '-'
+        else:
+            d = maxS
+            d_location = loc_maxS
+            d_sign = '+'
+
     g = gcd(n1, n2)
     n1g = n1 // g
     n2g = n2 // g
@@ -8098,7 +8167,7 @@ def ks_2samp(data1, data2, alternative='two-sided', method='auto'):
             prob = np.exp(expt)
 
     prob = np.clip(prob, 0, 1)
-    return KstestResult(d, prob)
+    return KstestResult(d, prob, statistic_location=d_location, statistic_sign=d_sign)
 
 
 def _parse_kstest_args(data1, data2, args, N):
@@ -8177,10 +8246,18 @@ def kstest(rvs, cdf, args=(), N=20, alternative='two-sided', method='auto'):
 
     Returns
     -------
-    statistic : float
-        KS test statistic, either D, D+ or D-.
-    pvalue : float
-        One-tailed or two-tailed p-value.
+    res: KstestResult
+        An object containing attributes:
+
+        statistic : float
+            KS test statistic, either D, D+ or D-.
+        pvalue : float
+            One-tailed or two-tailed p-value.
+
+        statistic_location: float
+            Value of the stochastic variable closest to the KS statistic. By
+            definition, the distance between the two sampled distributions at
+            this value is largest.
 
     See Also
     --------


Code review comments:
diff_hunk: @@ -8098,7 +8167,7 @@ def ks_2samp(data1, data2, alternative='two-sided', method='auto'):
             prob = np.exp(expt)
 
     prob = np.clip(prob, 0, 1)
-    return KstestResult(d, prob)
+    return KstestResult(d, prob, statistic_location=d_location, statistic_sign=d_sign)
path: scipy/stats/_stats_py.py
text: @user1:
```suggestion
    return KstestResult(d, prob, statistic_location=d_location,
                        statistic_sign=d_sign)
```
PEP8: `scipy/stats/_stats_py.py:8170:80: E501 line too long (86 > 79 characters)`

diff_hunk: @@ -7549,7 +7549,8 @@ def chisquare(f_obs, f_exp=None, ddof=0, axis=0):
                             lambda_="pearson")
 
 
-KstestResult = namedtuple('KstestResult', ('statistic', 'pvalue'))
+KstestResult = _make_tuple_bunch(['statistic', 'pvalue'],
path: scipy/stats/_stats_py.py
text: @user1:
```suggestion
KstestResult = _make_tuple_bunch('KstestResult', ['statistic', 'pvalue'],
```
I think this is causing most of the test failures, if not all.