How to write program that takes a vector and a CDF and performs the Kolmogorov-Smirnov test?

Question

I've been reading the Wikipedia articles about various statistical tests, and trying to code up programs that perform them. However, the article on the Kolmogorov-Smirnov test is rather unclear to me. Can somebody explain, in simple language, how I would write a program that takes a list of real numbers and a CDF, and computes a p-value for how closely they match?

Answer 1

Here is an answer to a related question that as part of the answer gives an example of computing the statistic and a graph demonstrating the general idea.

Answer 2

I can't make any sense of the code that ttnphns posted. But the linked algorithms document has helped.

First, it appears that if you take your $n$ samples and sort them into ascending order, then the emperical cumulative distribution function can be defined as

$$ \hat{F}(x) = \left\{ \begin{array}{ll} 0 & \text{if } x < x_0 \\ i/n & \text{if } x_i \leq x < x_{i+1} \\ 1 & \text{if } x_n < x \\ \end{array} \right. $$

(Assuming that the $x_i$ are your samples in ascending order.)

It now remains only to discover the maximum distance between $\hat{F}(x)$ and your target CDF $F(x)$.

It appears that finding the minimum distance would be quite hard, since it could be located anywhere in a given span. (The CDF might cross the EDF at any point.) However, since a CDF is, by definition, monotonically increasing, it follows that the maximum distance has to be at one end or the other of each span.

Edit: Corrected formulas. (Oops!)

So, it seems the complete algorithm is this:

• Take your samples $s_i$ and sort them into ascending order $x_i$.

• For $i=0$ to $n$, compute $d_L = |i/n - F(x_i)|$ and $d_H = |i/n - F(x_{i-1})|$.

• The largest $d_L$ or $d_H$ found is the KS statistic.

To get a p-value, you need to compare the statistic against the CDF of the Kolgomorov distribution. Apparently this involves the elliptic theta function or something... o_O

PS. It seems to not matter whether you choose $d_H = |i/n - F(x_{i-1})|$ or $d_H = |(i-1)/n - F(x_i)|$. Both do essentially the same thing.

Answer 3

Below I've written for you SPSS code of one-sample K-S test following SPSS algorithms document. In my example, an observed variable X is tested against uniform theoretical distribution. Despite that the snippet is not pseudocode you will understand easily what it is doing if you compare the comments with the document.

*Compute empirical cumulative distribution cum_e, and also sample size n.
rank X /rfraction into cum_e /n into n /ties= high /print= no .

*Let's obtain parameters for theoretical uniform distribution from the data (but one could specify arbitrary values).
aggregate /outfile= * mode= addvari /p1= min(X) /p2= max(X).

*Compute theoretical cumulative distribution.
compute cum_t= (X-p1)/(p2-p1).

*The two differences between the two cumulative curves.
*d1 are the differences when the empirical curve holds the right of the theoretical curve.
*d2 are the differences when the empirical curve holds the left of the theoretical curve.
*(they are needed both because theoretical distribution is continuous whereas empirical one is willy-nilly descrete).
*The utmost negative difference ("D-") is the minimum in d1.
*The utmost positive difference ("D+") is the maximum in d2.
compute casenum= $casenum.
sort cases by X.
compute d1= lag(cum_e)-cum_t.
compute d2= cum_e-cum_t.
sort cases by casenum.
aggregate /outfile= * mode= addvari over= yes /break= break /d1= min(d1) /d2= max(d2).
         /*Values "D-" (as d1 now) and "D+" (as d2 now)

*The test value, z.
do if $casenum=1.
-compute z= sqrt(n)*max(abs(d1),abs(d2)).

*Its 2-tailed significance sig.
-do if z>=0 and z<.27.
- compute sig= 1.
-else if z>=.27 and z<1.
- compute sig= exp(-1.233701*z**-2).
- compute sig= 1-(2.506628/z)*(sig+sig**9+sig**25).
-else if z>=1 and z<3.1.
- compute sig= exp(-2*z**2).
- compute sig= 2*(sig-sig**4+sig**9-sig**16).
-else if z>=3.1.
- compute sig= 0.
-end if.
end if.
execute.

Answer 4

Purely for giggles, here is an implementation of the KS test in Haskell:

module KS where

import Data.List (sort)

computeKS :: (Ord x, Num x) => (x -> Double) -> [x] -> Double
computeKS cdf us =
  let
    n  = fromIntegral (length us)
    xs = sort us
    ds = [ abs (i/n - cdf xi) `max` abs ((i-1)/n - cdf xi) | (i, xi) <- zip [0..] xs ]
  in maximum ds

Argue among yourselves whether this is any more readable...