Vectorizing Two Sample Kolmogorov-Smirnov test I'm implementing the Two-sample Kolmogorov-Smirnov test in MATLAB. I must admit I know very little of formal statistics and was simply trying to implement the description in wikipedia's page.

Right now I take as input the two vectors to be compared x1 and x2, and optionally the p-value that the check will run on. It outputs 1 if the vectors pass the test and 0 if they do not.

A couple of questions I have:


  Is there a better way to compute the empirical distribution functions? 


Right now I evaluate the empirical distribution function at the points defined in the vector t. In my naive approximation, I made it go from the smallest value of the vectors to the largest, using twice as many points as the vector with the most points. I have no idea is this is the best way to compute it, and suspect the vector doesn't have to be nearly as dense.


  Can I vectorize the search through the t vector?


Right now I have a for loop that goes through every value of tand computes the empirical distribution function in that way. It seems that there might be a way to vectorize this operation further and get rid of the for loop. Haven't figured out a good way to do it however.

And of course, any other suggestions are welcome!

My code is below:

function [OUT] = k_stest(x1,x2,p)

% Set default p-value
    if nargin == 2 
        p = 0.05;
    end 

% Compute lenghts
    N1 = length(x1); 
    N2 = length(x2);

% Set vector t over which the empirical distribution function will be computed.
    t = linspace(min(min([x1 x2])),max(max([x1 x2])),2*max([N1 N2]));

% Initialize the statistic, negative values guarantees it will be overwritten by first call.
    D = -1; 

% Set c from tabulated values
    if p == 0.10
        c = 1.22;
    elseif p == 0.05
        c = 1.36;
    elseif p == 0.025
        c = 1.48;
    elseif p == 0.01
        c = 1.63;
    elseif p == 0.005
        c = 1.73;
    elseif p == 0.001
        c = 1.95;
    else
        disp('Invalid p-value. Only p = 0.10 0.05 0.025 0.01 0.005 0.001 are supported')
        return
    end

% Search though the vector t, computing the empirical distribution function for each t, 
% and overwriting the statistic if a higher value is found.
    for i = 1:length(t)
        F1 = sum(x1&lt;=t(i))/N1;
        F2 = sum(x2&lt;=t(i))/N2;
        if abs(F2-F1) &gt; D
            D = abs(F2-F1);
        end
    end

% Compare the statistic to determine if the samples pass or fail.
    if D == -1;
       disp('Error, invalid input vectors')
       return
    elseif D &gt; c*sqrt((N1+N2)/(N1*N2));
        OUT = 0;
    else
        OUT = 1;
    end

end
I'll take this in chronological order, from the top down. I think I have covered all aspects of the code =)

function [OUT] = k_stest(x1,x2,p)


This looks nice, but it would be better if you had spaces between the input arguments. Also, k_stest is not a very good name, as it's hard to understand what it does just be looking at the name. I read in another review that "You don't need the brackets around the output variable, but you should always include it". I don't necessarily agree,  I usually omit them, but that's just a matter of personal preferences.



% Set default p-value
    if nargin == 2 
        p = 0.05;
    end 


This is OK. Depending on how you call the function you might add some additional checks. Is p a scalar? Are x1 and x2 vectors? This is not necessary, but can be useful if you want to use this function a long time from now and don't remember the specifics. assert can be very useful here.



% Compute lenghts
    N1 = length(x1); 
    N2 = length(x2);


numel (number of elements) is better than length. It's more robust and a lot faster for long vectors.



% Set vector t over which the empirical distribution function will be computed.
    t = linspace(min(min([x1 x2])),max(max([x1 x2])),2*max([N1 N2]));


This is dangerous, and a bit messy.

[x1 x2] assumes both x1 and x2 are row-vectors, i.e. horizontal vectors. A lot of functions in Matlab returns vectors as column-vectors i.e. vertical vectors by default, including the colon operator (:). 

If you happen to call this function with two vertical vectors instead of two horizontal vectors, this would fail (due to [x1 x2]).

Instead, you should do:

x12 = [x1(:); x2(:)]; 
t = linspace(min(x12), max(x12), 2*max([N1, N2]));


This way you're converting x1 and x2 to vertical vectors, and concatenate them vertically, prior to the linspace call. Now you only need one call to max and min.



The way you're initializing D is OK!



% Set c from tabulated values
    if p == 0.10
        c = 1.22;
    elseif p == 0.05
        c = 1.36;
    ...


Comparing floating point values using == is dangerous. For instance 0.1 + 0.2 == 0.3 will return false. Instead, you should compare the values with the desired value, but with some tolerance (for instance eps).

    if abs(p - 0.10) &lt; eps
        c = 1.22; 
    elseif abs(p-0.05) &lt; eps
        c = 1.36;


If two values are closer to each other than +/- eps then they are essentially equal.



% Search though the vector t, computing the empirical distribution function for each t, 
% and overwriting the statistic if a higher value is found.
    for i = 1:length(t)
        F1 = sum(x1&lt;=t(i))/N1;
        F2 = sum(x2&lt;=t(i))/N2;
        if abs(F2-F1) &gt; D
            D = abs(F2-F1);
        end
    end


You are right, vectorization is the way to go! Use bsxfun to create a logical matrix where each element that satisfy the condition x1 &lt;= t(ii) are true and the others are false. Then sum the entire mask, and take the maximum difference.

mask1 = bsxfun(@le, x1(:), t);
mask2 = bsxfun(@le, x2(:), t);
s1 = sum(mask1) / N1;
s2 = sum(mask2) / N2;
D = max(abs(s1 - s2));




% Compare the statistic to determine if the samples pass or fail.
    if D == -1;
       disp('Error, invalid input vectors')
       return
    elseif D &gt; c*sqrt((N1+N2)/(N1*N2));
        OUT = 0;
    else
        OUT = 1;
    end    
end


The check for invalid input vectors should be in the start of the function. I can't see a reason why you want it in the end instead. And as far as I can tell it will never fail (if the input is invalid then you will get an error further up).

Instead of the if and elseif, you can simply do:

OUT = D &lt;= c*sqrt((N1+N2)/(N1*N2));




To sum it all up, your function can be written like this. I've omitted the comment for simplicity. You should include them in your code.

function OUT = k_stest(x1, x2, p)

if nargin == 2
    p = 0.05;
end

% Add check for invalid input vectors!

N1 = numel(x1);
N2 = numel(x2);

x12 = [x1(:); x2(:)]; 
t = linspace(min(x12), max(x12), 2*max([N1, N2]));

if abs(p - 0.10) &lt; eps
    c = 1.22; 
elseif abs(p-0.05) &lt; eps
    c = 1.36;
% Continue with the rest of the values
%
end

mask1 = bsxfun(@le, x1(:), t);
mask2 = bsxfun(@le, x2(:), t);
s1 = sum(mask1) / N1;
s2 = sum(mask2) / N2;
D = max(abs(s1 - s2));

OUT = D &lt;= c*sqrt((N1+N2)/(N1*N2));

end