plotcdfkuiper(x,a,b,cdf,varargin)

function plotcdfkuiper(x,a,b,cdf,varargin)
% plot cumulative distribution for variates x between limits a,b and
% compare with the theoretical cdf using the Kuiper test.
%
% The Kuiper test is similar to the Kolmogorov-Smirnov test, but the
% Kolmogorov-Smirnov test is most sensitive around the median and less
% sensitive at the tails of the distribution, whereas the Kuiper test is
% equally sensitive to all values of x.
% Reference: WH Press, SA Teukolsky, WT Vetterling & BP Flannery, Numerical
% Recipes in Fortran, 2nd edition, Cambridge University Press 1992, Ch 14.3
%
% cdf = handle for cumulative distribution function
% varargin contains parameters for cdf
% plot limits and cdf are optional inputs
%
% example using normal distribution:
% cdf = @(x) 0.5*erfc(-x/sqrt(2));
% plotcdfkuiper(randn(1,100),-3,3,cdf)
%
% S Bocquet DSTO 7 August 2008 Tested with MATLAB Version 7.6.0.324 (R2008a)
%
n = numel(x);
x = reshape(x,1,n); % ensure x is a row vector
yn = (0:n)/n;
y2 = reshape(repmat(yn,2,1),1,2*n+2); % cumulative probability
x2 = reshape(repmat(sort(x),2,1),1,2*n); % sort and duplicate variates
% make sure plot limits are sensible
if nargin < 2 || isempty(a) || ~isfinite(a) || a > x2(end); a = floor(x2(1)); end;
if nargin < 3 || isempty(b) || ~isfinite(b) || b < a; b = ceil(x2(end)); end;
x2 = [-Inf x2 Inf];
plot(x2,y2,'b')
axis([a b 0 1])
title('Cumulative distribution')
xlabel('x')
ylabel('F(x)')
if nargin > 3
  % calculate Kuiper statistic and significance
  yc = cdf(x2(2:end-1),varargin{:});
  Kstat = max(y2(2:end-1)-yc) + max(yc-y2(2:end-1));
  prob = kuiperprob((sqrt(n)+0.155+0.24/sqrt(n)).*Kstat);
  label = char(['Kuiper statistic: ' num2str(Kstat,'%8.3G')], ...
    ['Significance:     ' num2str(prob,'%8.3G')]);
  xs = a:(b-a)/1000:b;
  hold on
  plot(xs,cdf(xs,varargin{:}),'k') % plot theoretical cdf
  h = text(0,0,label,'EdgeColor','k');
  ext = get(h,'Extent');
  set(h,'Position',[b-ext(3)-0.02*(b-a),0.1]); % put label in SE corner
  hold off
end
end

function QK = kuiperprob(lambda)
% probability for Kuiper statistic
if lambda < 0.3
  QK = 1;
elseif lambda > 5
  QK = 0;
else
  j = 1:15; % use ceil(sqrt(-log(eps)/2)/0.3)=15 terms to evaluate series
  ej = 2*(lambda*j).^2;
  QK = 2*sum((2*ej-1).*exp(-ej));
end
end