Home > atmlab > circular > circ_wwtest.m

circ_wwtest

PURPOSE ^

[pval, table] = circ_wwtest(alpha, idx, [w])

SYNOPSIS ^

function [pval table] = circ_wwtest(varargin)

DESCRIPTION ^

 [pval, table] = circ_wwtest(alpha, idx, [w])
 [pval, table] = circ_wwtest(alpha1, alpha2, [w1, w2])
   Parametric Watson-Williams multi-sample test for equal means. Can be
   used as a one-way ANOVA test for circular data.

   H0: the s populations have equal means
   HA: the s populations have unequal means

   Note: 
   Use with binned data is only advisable if binning is finer than 10 deg.
   In this case, alpha is assumed to correspond
   to bin centers.

   The Watson-Williams two-sample test assumes underlying von-Mises
   distributrions. All groups are assumed to have a common concentration
   parameter k.

   Input:
     alpha   angles in radians
     idx     indicates which population the respective angle in alpha
             comes from, 1:s
     [w      number of incidences in case of binned angle data]

   Output:
     pval    p-value of the Watson-Williams multi-sample test. Discard H0 if
             pval is small.
     table   cell array containg the ANOVA table

 PHB 3/19/2009

 References:
   Biostatistical Analysis, J. H. Zar

 Circular Statistics Toolbox for Matlab

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

circ_wwtest.m

SOURCE CODE ^

0001 function [pval table] = circ_wwtest(varargin)
0002 % [pval, table] = circ_wwtest(alpha, idx, [w])
0003 % [pval, table] = circ_wwtest(alpha1, alpha2, [w1, w2])
0004 %   Parametric Watson-Williams multi-sample test for equal means. Can be
0005 %   used as a one-way ANOVA test for circular data.
0006 %
0007 %   H0: the s populations have equal means
0008 %   HA: the s populations have unequal means
0009 %
0010 %   Note:
0011 %   Use with binned data is only advisable if binning is finer than 10 deg.
0012 %   In this case, alpha is assumed to correspond
0013 %   to bin centers.
0014 %
0015 %   The Watson-Williams two-sample test assumes underlying von-Mises
0016 %   distributrions. All groups are assumed to have a common concentration
0017 %   parameter k.
0018 %
0019 %   Input:
0020 %     alpha   angles in radians
0021 %     idx     indicates which population the respective angle in alpha
0022 %             comes from, 1:s
0023 %     [w      number of incidences in case of binned angle data]
0024 %
0025 %   Output:
0026 %     pval    p-value of the Watson-Williams multi-sample test. Discard H0 if
0027 %             pval is small.
0028 %     table   cell array containg the ANOVA table
0029 %
0030 % PHB 3/19/2009
0031 %
0032 % References:
0033 %   Biostatistical Analysis, J. H. Zar
0034 %
0035 % Circular Statistics Toolbox for Matlab
0036 
0037 % By Philipp Berens, 2009
0038 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
0039 
0040 [alpha, idx, w] = processInput(varargin{:});
0041 
0042 % number of groups
0043 u = unique(idx);
0044 s = length(u);
0045 
0046 % number of samples
0047 n = sum(w);
0048 
0049 % compute relevant quantitites
0050 pn = zeros(s,1); pr = pn;
0051 for t=1:s
0052   pidx = idx == u(t);
0053   pn(t) = sum(pidx.*w);
0054   pr(t) = circ_r(alpha(pidx),w(pidx));
0055 end
0056 
0057 r = circ_r(alpha,w);
0058 rw = sum(pn.*pr)/n;
0059 
0060 % make sure assumptions are satisfied
0061 checkAssumption(rw,mean(pn))
0062 
0063 % test statistic
0064 kk = circ_kappa(rw);
0065 beta = 1+3/(8*kk);    % correction factor
0066 A = sum(pr.*pn) - r*n;
0067 B = n - sum(pr.*pn);
0068 
0069 F = beta * (n-s) * A / (s-1) / B;
0070 pval = 1 - fcdf(F,s-1,n-s);
0071 
0072 na = nargout;
0073 if na < 2
0074   printTable;
0075 end
0076 prepareOutput;
0077 
0078 
0079   function printTable
0080     
0081     fprintf('\nANALYSIS OF VARIANCE TABLE (WATSON-WILLIAMS TEST)\n\n');
0082     fprintf('%s\t\t\t\t%s\t%s\t\t%s\t\t%s\t\t\t%s\n', ' ' ,'d.f.', 'SS', 'MS', 'F', 'P-Value');
0083     fprintf('--------------------------------------------------------------------\n');
0084     fprintf('%s\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', 'Columns', s-1 , A, A/(s-1), F, pval);
0085     fprintf('%s\t\t%u\t\t%.2f\t%.2f\n', 'Residual ', n-s, B, B/(n-s));
0086     fprintf('--------------------------------------------------------------------\n');
0087     fprintf('%s\t\t%u\t\t%.2f', 'Total   ',n-1,A+B);
0088     fprintf('\n\n')
0089 
0090   end
0091 
0092   function prepareOutput
0093     
0094     if na > 1
0095       table = {'Source','d.f.','SS','MS','F','P-Value'; ...
0096                'Columns', s-1 , A, A/(s-1), F, pval; ...
0097                'Residual ', n-s, B, B/(n-s), [], []; ...
0098                'Total',n-1,A+B,[],[],[]};
0099     end
0100   end
0101 end
0102 
0103 
0104 
0105 function checkAssumption(rw,n)
0106 
0107   if n > 10 && rw<.45
0108     warning('Test not applicable. Average resultant vector length < 0.45.') %#ok<WNTAG>
0109   elseif n > 6 && rw<.5
0110     warning('Test not applicable. Average number of samples per population < 11 and average resultant vector length < 0.5.') %#ok<WNTAG>
0111   elseif n >=5 && rw<.55
0112     warning('Test not applicable. Average number of samples per population < 7 and average resultant vector length < 0.55.') %#ok<WNTAG>
0113   elseif n < 5
0114     warning('Test not applicable. Average number of samples per population < 5.') %#ok<WNTAG>
0115   end
0116 
0117 end
0118 
0119 
0120 function [alpha, idx, w] = processInput(varargin)
0121 
0122   if nargin == 4
0123     alpha1 = varargin{1}(:);
0124     alpha2 = varargin{2}(:);
0125     w1 = varargin{3}(:);
0126     w2 = varargin{4}(:);
0127     alpha = [alpha1; alpha2];
0128     idx = [ones(size(alpha1)); ones(size(alpha2))];
0129     w = [w1; w2];
0130   elseif nargin==2 && sum(abs(round(varargin{2})-varargin{2}))>1e-5
0131     alpha1 = varargin{1}(:);
0132     alpha2 = varargin{2}(:);
0133     alpha = [alpha1; alpha2];
0134     idx = [ones(size(alpha1)); 2*ones(size(alpha2))];
0135     w = ones(size(alpha));
0136   elseif nargin==2
0137     alpha = varargin{1}(:);
0138     idx = varargin{2}(:);
0139     if ~(size(idx,1)==size(alpha,1))
0140       error('Input dimensions do not match.')
0141     end
0142     w = ones(size(alpha));  
0143   elseif nargin==3
0144     alpha = varargin{1}(:);
0145     idx = varargin{2}(:);
0146     w = varargin{3}(:);
0147     if ~(size(idx,1)==size(alpha,1))
0148       error('Input dimensions do not match.')
0149     end
0150     if ~(size(w,1)==size(alpha,1))
0151       error('Input dimensions do not match.')
0152     end  
0153   else
0154     error('Invalid use of circ_wwtest. Type help circ_wwtest.')
0155   end
0156 end
0157

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