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
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
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<=t(i))/N1; F2 = sum(x2<=t(i))/N2; if abs(F2-F1) > 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 > c*sqrt((N1+N2)/(N1*N2)); OUT = 0; else OUT = 1; end end