+
−% KSpValue Compute p-Value of the Kolmogorov - Smirnov distribution
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−%
+
−% Compute p-Value of the Kolmogorov - Smirnov distribution
+
−%
+
−% CALL
+
−%
+
−% pValue = utils.math.KSpValue(KSstatistic,n1,n2);
+
−%
+
−%
+
−% INPUT
+
−%
+
−% - KSstatistic, value of the statistic of the KS distribution.
+
−% Corresponding at KSstatistic = max(abs(CD1-CD2))
+
−% - length of the first data series
+
−% - length of the second data series
+
−%
+
−% References:
+
−%   Marsaglia, G., W.W. Tsang, and J. Wang (2003), "Evaluating Kolmogorov`s
+
−%         Distribution", Journal of Statistical Software, vol. 8, issue 18.
+
−%
+
−%
+
−% L Ferraioli 06-12-2010
+
−%
+
−% $Id: KSpValue.m,v 1.3 2011/07/14 07:10:16 mauro Exp $
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−function pValue = KSpValue(KSstatistic, n1, n2)
+
−  
+
−  if isempty(n2)
+
−    n = n1; % test against theoretical distribution
+
−  else
+
−    n = n1*n2/(n1+n2); % test of two empirical distributions
+
−  end
+
−  s = n*KSstatistic^2;
+
−  
+
−  % Following the recipe described in described in Marsaglia, et al.
+
−  % For d values that are in the far tail of the distribution (i.e.
+
−  % p-values > .999), the following lines will speed up the computation
+
−  % significantly, and provide accuracy up to 7 digits.
+
−  if s == 0
+
−    pValue = 0;
+
−  else
+
−    if (s > 7.24) || ((s > 3.76) && (n > 99))
+
−      pValue = 2*exp(-(2.000071+.331/sqrt(n)+1.409/n)*s);
+
−    else
+
−      % Express d as d = (k-h)/n, where k is a +ve integer and 0 < h < 1.
+
−      k = ceil(KSstatistic*n);
+
−      h = k - KSstatistic*n;
+
−      m = 2*k-1;
+
−      
+
−      % Create the H matrix, which describes the CDF, as described in Marsaglia,
+
−      % et al.
+
−      if m > 1
+
−        c = 1./gamma((1:m)' + 1);
+
−        
+
−        r = zeros(1,m);
+
−        r(1) = 1;
+
−        r(2) = 1;
+
−        
+
−        T = toeplitz(c,r);
+
−        
+
−        T(:,1) = T(:,1) - (h.^[1:m]')./gamma((1:m)' + 1);
+
−        
+
−        T(m,:) = fliplr(T(:,1)');
+
−        T(m,1) = (1 - 2*h^m + max(0,2*h-1)^m)/gamma(m+1);
+
−      else
+
−        T = (1 - 2*h^m + max(0,2*h-1)^m)/gamma(m+1);
+
−      end
+
−      
+
−      % Scaling before raising the matrix to a power
+
−      if ~isscalar(T)
+
−        lmax = max(eig(T));
+
−        T = (T./lmax)^n;
+
−      else
+
−        lmax = 1;
+
−      end
+
−      
+
−      % Pr(Dn < d) = n!/n * tkk ,  where tkk is the kth element of Tn = T^n.
+
−      % p-value = Pr(Dn > d) = 1-Pr(Dn < d)
+
−      pValue = (1 - exp(gammaln(n+1) + n*log(lmax) - n*log(n)) * T(k,k));
+
−    end
+
−  end
+
−  pValue = abs(pValue);
+
−  
+
−end