# what is the algorithm of getting maximum $D$ value in ks.test?

When I was writing R code to plot two empirical cumulative density curves, I came across ks.test in R. I looked into the code of ks.test and ran into this question.

I am wondering whether anyone can explain to me why this expression

z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))


gives the maximum $D$ values, specifically why order(w)<=n.x? Thanks!

The order function tells you the index position of the smallest, second smallest, and so on up to largest. So if order(w)[1] is 8, the eighth value in w is the smallest.

The vector $w$ is the combined sample (it has all the $x$ values, then all the $y$ values). Since all the $x$ values are first, when the position is $\leq$ the $x$-sample size (order(w)<=n.x) we're looking at that $w$ value being from the $x$-sample.

So when order(w) is less than or equal to the number of values in $x$, the ecdf of $x$ undergoes a jump at that value of $w$, otherwise the ecdf of $y$ does. By taking negative jumps for the $y$ and positive for the $x$, we are looking at the changes in the difference between the ecdfs. Adding those jumps up (via cumsum) then gives the difference in ecdfs, which is what that line computes.

You could do the more obvious z <- ecdf(x)(sort(w))-ecdf(y)(sort(w)) instead, but it would be slower.

This should help. Basically, KS test compares cumulative probabilities. In order to construct the CDF, you need to order the observations. Then you go from the smallest to largest, and keep counting. When you reach the largest observation your count is 100%, i.e. $P=\frac{i}{N}$ where $N$ is the number of observations, and $i$ is the index of this observation. While doing so you compare $P$ with the theoretical CDF, and get the maximum deviation

• +1 You might want to expand your answer to the two sample case. Dec 11, 2014 at 0:07
• The two sample case is very similar, but instead of theoretical CDF, the second sample's empirical CDF is used. I think that the test statistics is slightly different though. Dec 11, 2014 at 4:10

This is clever code; it deserves to be understood well. I will therefore unravel it carefully, explaining each part. As with much functional code, we have to work from the inside out.

The calculation is apropros the two-sample Kolmogorov-Smirnov test. One sample has been stored in array $\mathrm{x}=(x_1,x_2,\ldots, x_n)$ and the other in $\mathrm{y}=(y_1,y_2,\ldots,y_m)$ (where to avoid towers of subscripts I write $n$ for n.x and $m$ for n.y). The array $\mathrm w$ is their concatenation

$$\mathrm{w} = (x_1, x_2, \ldots, x_n,\, y_1, y_2, \ldots, y_m).$$

The purpose of order is to give the indexes of $\mathrm{w}$ in increasing sequence. Writing the output of order as $\sigma = (\sigma(1), \sigma(2), \ldots, \sigma(n+m))$, this means that

$$w_{\sigma(i)} \le w_{\sigma(j)}$$

whenever $i \lt j$.

The expression order(w) <= n.x creates an indicator vector with only zeros and ones for its entries. It will have ones exactly at those places where $\sigma(i)\le n$. Since $\sigma(i)$ is an entry into $\mathrm w$, and the first $n$ entries of $\mathrm w$ come from $\mathrm x$, this is tantamount to creating an indicator of where the elements of $\mathrm x$ would be after sorting $\mathrm w$.

The effect of cumsum(ifelse is to use this indicator to walk around on a vertical line starting at zero. Every time an element of $\mathrm x$ is encountered, move $1/n$ up; every time an element of $\mathrm y$ is encountered, move $1/m$ down.

Compare this recipe to the following illustration of the empirical distribution functions (posted in an answer to a related question at https://stackoverflow.com/questions/27233738, but with the roles of $\mathrm x$ and $\mathrm y$ reversed there):

The red graph is the empirical CDF (ECDF) of $\mathrm{x}$: by construction, as you scan from left to right, it leaps upward by $1/n$ each time an element of $\mathrm{x}$ is encountered on the horizontal axis. Similarly, the blue graph is the empirical CDF of $\mathrm{y}$; its leaps occur in quantum units of $1/m$. The vertical black line shows one stage of this left-to-right scan; its vertical length represents the difference between the red graph and the blue graph. (If the red graph drops below the blue graph, the length is considered negative.)

We can know the length of the vertical line in a simple fashion. At the outset, it has zero length. Every time we scan across an element of $\mathrm {x}$ it lengthens by $1/n$, but every time we scan across an element of $\mathrm{y}$ it is shortened by $1/m$. As we saw above, that's precisely what the magic R code does!

Consequently,

cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))

is an array that contains all the distinct (signed) lengths attained by the vertical black line throughout the scan. Assuming there are no ties between some value of $\mathrm x$ and some value of $\mathrm y$ (which is an assumption of this test which is explicitly checked in the code beforehand), its maximum is the greatest distance reached by the ECDF of $\mathrm{x}$ above the ECDF of $\mathrm {y}$, while its minimum is the negative of the greatest distance reached by the ECDF of $\mathrm{x}$ below the ECDF of $\mathrm {y}$. Their absolute values are precisely the two Kolmogorov-Smirnov statistics $D^{+}$ and $D^{-}$ and the larger of those two absolute values is the K-S statistic $D$.

• +1!!! Question, though: you write "every time an element of y is encountered, move 1/m down." Should that be "...move 1/m up?" In the graph below your quoted statement scanning both x and y from left to right, both eCDFs are monotonically increasing (x by 1/n, y by 1/m). Or am I misunderstanding something? (Am sure I am, just puzzling.) Dec 11, 2014 at 6:26
• Trying to answer my own question reading a little farther... by "walk around on a vertical line" you mean something like "track the vertical distance between x's eCDF and y's eCDF?" Dec 11, 2014 at 6:32
• @Alexis Yes: that is how you should interpret the location on the vertical line.
– whuber
Dec 11, 2014 at 14:30