Home > atmlab > circular > circ_kuipertest.m

circ_kuipertest

PURPOSE ^

[pval, k, K] = circ_kuipertest(sample1, sample2, res, vis_on)

SYNOPSIS ^

function [pval, k, K] = circ_kuipertest(alpha1, alpha2, res, vis_on)

DESCRIPTION ^

 [pval, k, K] = circ_kuipertest(sample1, sample2, res, vis_on)

   The Kuiper two-sample test tests whether the two samples differ 
   significantly.The difference can be in any property, such as mean 
   location and dispersion. It is a circular analogue of the 
   Kolmogorov-Smirnov test.  
 
   H0: The two distributions are identical.
   HA: The two distributions are different.

 Input: 
   alpha1    fist sample (in radians)
   alpha2    second sample (in radians)
   res       resolution at which the cdf is evaluated
   vis_on    display graph

 Output:
   pval        p-value; the smallest of .10, .05, .02, .01, .005, .002,
               .001, for which the test statistic is still higher
               than the respective critical value. this is due to
               the use of tabulated values. if p>.1, pval is set to 1.
   k           test statistic
   K           critical value
 
 References:
   Batschelet, 1980, p. 112

 Circular Statistics Toolbox for Matlab

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

circ_kuipertest.m

SOURCE CODE ^

0001 function [pval, k, K] = circ_kuipertest(alpha1, alpha2, res, vis_on)
0002 
0003 % [pval, k, K] = circ_kuipertest(sample1, sample2, res, vis_on)
0004 %
0005 %   The Kuiper two-sample test tests whether the two samples differ
0006 %   significantly.The difference can be in any property, such as mean
0007 %   location and dispersion. It is a circular analogue of the
0008 %   Kolmogorov-Smirnov test.
0009 %
0010 %   H0: The two distributions are identical.
0011 %   HA: The two distributions are different.
0012 %
0013 % Input:
0014 %   alpha1    fist sample (in radians)
0015 %   alpha2    second sample (in radians)
0016 %   res       resolution at which the cdf is evaluated
0017 %   vis_on    display graph
0018 %
0019 % Output:
0020 %   pval        p-value; the smallest of .10, .05, .02, .01, .005, .002,
0021 %               .001, for which the test statistic is still higher
0022 %               than the respective critical value. this is due to
0023 %               the use of tabulated values. if p>.1, pval is set to 1.
0024 %   k           test statistic
0025 %   K           critical value
0026 %
0027 % References:
0028 %   Batschelet, 1980, p. 112
0029 %
0030 % Circular Statistics Toolbox for Matlab
0031 
0032 % By Marc J. Velasco and Philipp Berens, 2009
0033 % velasco@ccs.fau.edu
0034 
0035 
0036 if nargin < 3
0037     res = 100;
0038 end
0039 if nargin < 4
0040     vis_on = 0;
0041 end
0042 
0043 n = length(alpha1(:));
0044 m = length(alpha2(:));
0045 
0046 % create cdfs of both samples
0047 [phis1 cdf1 phiplot1 cdfplot1] = circ_samplecdf(alpha1, res);
0048 [phis2 cdf2 phiplot2 cdfplot2] = circ_samplecdf(alpha2, res);
0049 
0050 % maximal difference between sample cdfs
0051 [dplus, gdpi] = max([0 cdf1-cdf2]);
0052 [dminus, gdmi] = max([0 cdf2-cdf1]);
0053 
0054 % calculate k-statistic
0055 k = n * m * (dplus + dminus);
0056 
0057 % find p-value
0058 [pval K] = kuiperlookup(min(n,m),k/sqrt(n*m*(n+m)));
0059 K = K * sqrt(n*m*(n+m));
0060 
0061 
0062 % visualize
0063 if vis_on
0064     figure 
0065     plot(phiplot1, cdfplot1, 'b', phiplot2, cdfplot2, 'r');
0066     hold on
0067     plot([phis1(gdpi-1), phis1(gdpi-1)], [cdf1(gdpi-1) cdf2(gdpi-1)], 'o:g');
0068     plot([phis1(gdmi-1), phis1(gdmi-1)], [cdf1(gdmi-1) cdf2(gdmi-1)], 'o:g');
0069     hold off
0070     set(gca, 'XLim', [0, 2*pi]);
0071     set(gca, 'YLim', [0, 1.1]);
0072     xlabel('Circular Location')
0073     ylabel('Sample CDF')
0074     title('CircStat: Kuiper test')
0075     h = legend('Sample 1', 'Sample 2', 'Location', 'Southeast');
0076     set(h,'box','off')
0077     set(gca, 'XTick', pi*0:.25:2)
0078     set(gca, 'XTickLabel', {'0', '', '', '', 'pi', '', '', '', '2pi'}) 
0079 end
0080 
0081 
0082 
0083 end
0084 
0085 function [p K] = kuiperlookup(n, k)
0086 
0087 load kuipertable.mat;
0088 alpha = [.10, .05, .02, .01, .005, .002, .001];
0089 nn = ktable(:,1); %#ok<COLND>
0090 
0091 % find correct row of the table
0092 [easy row] = ismember(n, nn);
0093 if ~easy
0094    % find closest value if no entry is present)
0095    row = length(nn) - sum(n<nn); 
0096    if row == 0
0097        error('N too small.');
0098    else
0099       warning('N=%d not found in table, using closest N=%d present.',n,nn(row)) %#ok<WNTAG>
0100    end
0101 end
0102 
0103 % find minimal p-value and test-statistic
0104 idx = find(ktable(row,2:end)<k,1,'last');
0105 if ~isempty(idx)
0106   p = alpha(idx);
0107 else
0108   p = 1;
0109 end
0110 K = ktable(row,idx+1);
0111 
0112 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005