2

I'm comparing two different empirical cumulative distribution functions using the KS-test, and I'd like to extract the location (in the ECDF) where the maximum of the test statistic is.

Question: Using R, is there a convenient way to extract that, perhaps from the ks.test function or otherwise?

Thanks for any and all comments.

5
  • Are you thinking the test statistic is a pointwise function? Seems to me that the stats.stackexchange.com moderators are shirking their responsibilities in failing to clarify the statistical confusion.
    – IRTFM
    Commented Dec 1, 2014 at 17:58
  • @Bonded On the contrary, this is a purely programming question with a purely programming answer. That makes it clearly off-topic on CV; whether it is considered on-topic here is a matter for this community to decide. To help with that decision I have posted an answer that consists only of code and includes no statistical analysis or reasoning.
    – whuber
    Commented Dec 1, 2014 at 18:08
  • Admittedly, I have found that the proper location for asking R questions (stats or stackoverflow) can be rather muddied. Whuber is correct in that my central confusion was in the coding rather than the statistics. Thank you, Whuber, as that does answer what I was looking for. Looking back, I was getting hung up on how to go from the max statistic to the index, or the function which.max().
    – Ashe
    Commented Dec 1, 2014 at 18:34
  • 1
    Questions to SO shouldn't just ask for code; they should present an example or test case in code and lay out a plan for constructing an answer.
    – IRTFM
    Commented Dec 1, 2014 at 19:26
  • @Bonded Fair enough--and we certainly do try to apply that criterion when making migration decisions over at CV. In what perhaps is a somewhat generous spirit, one might take the guidance in this question (such as the idea to exploit ks.test) as a plan for obtaining the answer (and indeed it turned out to be a good idea). It's difficult to conceive of how one would present any kind of useful code to represent partial progress in this case.
    – whuber
    Commented Dec 1, 2014 at 23:09

1 Answer 1

7

It does not appear you can extract such a location (which might not be unique, BTW) from the output of ks.test, but by emulating the key calculation there you can obtain the answer:

compare <- function(x, y) {
  n <- length(x); m <- length(y)
  w <- c(x, y)
  o <- order(w)
  z <- cumsum(ifelse(o <= n, m, -n))
  i <- which.max(abs(z))
  w[o[i]]
}

The calculation through z <- ... is from the ks.test source, while the last two lines (fairly clearly) find the location where the maximum deviation is attained.

As an example, let's generate two datasets and compare them:

set.seed(17)
x <- rnorm(30)
y <- rnorm(20, sd=2/3)
u <- compare(x,y)

The reported value of u is 0.04946235. To see whether this is correct, check it against the ECDFs and the output of ks.test:

e.x <- ecdf(x)
e.y <- ecdf(y)
abs(e.x(u) - e.y(u))
ks.test(x,y)$statistic

The output in both cases is 0.4166667, indicating perfect agreement. A plot of the situation will clarify what is going on:

plot(e.x, col="Blue", main="ECDF", xlab="Value", ylab="Probability")
plot(e.y, add=TRUE, col="Red")
lines(c(u,u), c(0,1), col="Gray")
lines(c(u,u), c(e.x(u), e.y(u)), lwd=2)

It shows both ECDFs and marks the location found by compare (namely, u) with a vertical line: it is supposed to indicate the place where the two graphs attain their greatest vertical separation.

Figure

6
  • 1
    I think your answer is an excellent R example, but I also think that the OP did not understand the statistical need for construction of a common (sorted) range c(x,y) as arguments for comparison of the two different ECDFs. I would have set the plot limits to xlim=range(c(x,y)) to more generally handle the problem of 2 sets with ranges that don't overlap as well.
    – IRTFM
    Commented Dec 1, 2014 at 18:49
  • @Bonded Thank you for the good suggestions. (I didn't worry about plot limits because the plot is purely to illustrate what is going on. Some of the ks.test source code is devoted to handling the case of non-overlap, BTW.) I cannot speak to what the OP might or might not have understood, but only hope that the code and the plot will help all readers appreciate the question.
    – whuber
    Commented Dec 1, 2014 at 19:14
  • I am quite a noobie but why does the ks.test doesn't the correction for large $n$? For large $n$ the code part form you needs to be changed to: (sqrt((n * m) / (n+m)) ) * abs(e.x(u) - e.y(u)) (imho). If It do this, the test statistic is different from the one in given by $ks.test(x,y)$statistic$. I habe $n$ which sizes around 600k. Imho that is a large $n$. Hence, there is a discrepancy I do not understand
    – kn1g
    Commented Jun 2, 2020 at 13:31
  • 1
    @kn1g First, you appear to be confusing an asymptotic approximation with the statistic. Second, the KS test is useless for datasets that size. Search Cross Validated for extensive discussions of this test and Normality testing generally.
    – whuber
    Commented Jun 2, 2020 at 13:37
  • Thanks. I try to compare to return distributions. I started googling when the KS test gets useless because I didn't know that. I understand the point with my confusion of asymptotic approximation and the statistic. Thanks for pointing this out. Also, I am still looking for a better way to compare the distributions as they are returns and therefore not i.i.d.
    – kn1g
    Commented Jun 2, 2020 at 13:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.