%
% normtest - test the cumulative distribution of a
%   given data set against the normal distribution
%   and apply the Kolmogorov-Smirnov Test
%
%   [ pks, h ] = normtest( data )
%
% where pks is probability data is normal (0 < pks < 1)
%
% Bob Newell, Sept 1997
%
function [ pks, h ] = normtest( data )
%
% parameters for normal cumulative prob approximation
%  from C.Hastings Jr, Princeton Univ Press, 1955
%
d = [ 0.049867347 0.0211410061 0.0032776263 ...
      0.0000380036 0.0000488906 0.000005383 ];
%
% calculate actual cumulative probability function
%
m = 100; if length(data) < m, m = length(data); end;
[ n, x ] = hist( (data-mean(data))/std(data), m );
f = cumsum(n) / sum(n);
%
% calculate normal cumulative probability function
%
xx = ones( size(x) );
ss = xx;
for i = 1:6,
  xx = xx .* abs(x);
  ss = ss + d(i) * xx;
end;
ksd = 0;
for i = 1:length(x),
  if x(i) > 0,
    px(i) = 1 - 0.5 * ss(i) ^ (-16);
  else
    px(i) = 0.5 * ss(i) ^ (-16);
  end;
% this is the K-S metric calculation
  d = abs( f(i) - px(i) );
  if d > ksd, ksd = d; end;
%
end;
%
fontsz = 16;
h = figure;
plot( px, f, '*' )
line( [ 0 1 ], [ 0 1 ] );
set( gca, 'fontsize', fontsz );
xlabel( 'Normal cumulative probability' );
ylabel( 'Actual cumulative probability' );
%
% calculate the K-S probability
%
a2 = - 2 * m * ksd * ksd;
fac = 2;
pks = 0;
for i = 1:100,
  pks = pks + fac * exp( a2 * i * i );
  fac = - fac;
end;
%
ht = text( 0.05, 0.90, 'Kolmogorov-Smirnov' );
set( ht, 'fontsize', fontsz );
ht = text( 0.10, 0.80, [ 'Probability = ' int2str(round(100*pks)) ' %' ] );
set( ht, 'fontsize', fontsz );
%