From octave-sources-request at bevo dot che dot wisc dot edu Thu Apr 10 07:51:14 2003
Subject: Calculation of the Smirnov-Kolmogorov table
From: Alexandre Vial <alexandre dot vial at utt dot fr>
To: octave-sources at bevo dot che dot wisc dot edu
CC: jrv at vanzandt dot mv dot com
Date: Thu, 10 Apr 2003 14:49:45 +0200

Almost two years ago, James R Van Zandt wrote:
(http://www.octave.org/mailing-lists/octave-sources/2001/100)
 >
 >> What is the relationship between the distribution calculated by
 >> Octave's function and the tabulated distribution?  I'd like to add
 >> that to the documentation.

Actually, what kolmogorov_smirnov_cdf(alpha) gives is 
P(sqrt(n).Dn<alpha), and this is a valid approximation for n>100 (some 
books say 35, but then the approximation is bad).
For small value of n, you have to use the following code (WARNING: it 
gives P(Dn>alpha), not P(Dn<alpha), for big values of n, you have to 
compare the results of this function with 
1-kolmogorov_smirnov_cdf(sqrt(n)*alpha) ) :

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% skbl.m  (Smirnov-Kolmogorov Bi-Lateral)
% this function needs the number of values "n", the max difference
% "alpha"
% it returns p_ks=P(Dn > alpha)
function p_ks=skbl(n,alpha)

   kmax=n*(1-alpha);
   k=0:floor(kmax);
   t1=gamma(n+1)./(gamma(k+1).*gamma(n-k+1));
   t2=(alpha+k/n).^(k-1);
   t3=(1-alpha-k/n).^(n-k);
   p_ks=sum(2*t1*alpha.*t2.*t3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Taking the example of James R Van Zandt (n=7, Dn=0.381) :
octave:2> skbl(7,0.381)
ans = 0.20114
we obtain 0.2, what was expected

Now, if you want to calculate the SK table (as given here for example 
http://www.eridlc.com/onlinetextbook/appendix/table7.htm with slightly 
different level of significance), you have to use a dichotomy algorithm 
(probably not optimized, but it works for me):

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p_line=[0.2 0.1 0.05 0.025 0.02 0.01 0.005];
n_max=100;
n_line=1:n_max;
tableau=zeros(n_max+1,8);

tableau(1,2:8)=p_line;

for n=n_line
	tableau(n+1,1)=n;
	colonne=1;
	for p=p_line
		colonne=colonne+1;
		dn=1; % max value
		delta=0.5;
		res=skbl(n,dn);
		while (res>1.0001*p  |  res<0.9999*p)
			if (res>1.0001*p)
				dn=dn+delta;
			end
			if (res<0.9999*p)
				dn=dn-delta;
			end
			delta=delta/2;
			res=skbl(n,dn);			
		end
		tableau(n+1,colonne)=dn;
	end
	n
	tableau(n+1,2:8)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I hope you will enjoy it !

Alexandre Vial

-- 
Alexandre Vial - Université de Technologie de Troyes
Laboratoire de Nanotechnologie et d'Instrumentation Optique
CNRS FRE 2671		         Tel : +33 (0)3 25 71 80 28
12 rue Marie Curie BP-2060	 Fax : +33 (0)3 25 71 56 75
F-10010 TROYES Cedex		 http://www-lnio.utt.fr