forked from circstat/circstat-matlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
circ_cmtest.m
90 lines (71 loc) · 2.08 KB
/
circ_cmtest.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
function [pval med P] = circ_cmtest(varargin)
%
% [pval, med, P] = circ_cmtest(alpha, idx)
% [pval, med, P] = circ_cmtest(alpha1, alpha2)
% Non parametric multi-sample test for equal medians. Similar to a
% Kruskal-Wallis test for linear data.
%
% H0: the s populations have equal medians
% HA: the s populations have unequal medians
%
% Input:
% alpha angles in radians
% idx indicates which population the respective angle in alpha
% comes from, 1:s
%
% Output:
% pval p-value of the common median multi-sample test. Discard H0 if
% pval is small.
% med best estimate of shared population median if H0 is not
% discarded at the 0.05 level and NaN otherwise.
% P test statistic of the common median test.
%
%
% PHB 7/19/2009
%
% References:
% Fisher NI, 1995
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% [email protected] - www.kyb.mpg.de/~berens/circStat.html
[alpha, idx] = processInput(varargin{:});
% number of groups
u = unique(idx);
s = length(u);
% number of samples
N = length(idx);
% total median
med = circ_median(alpha);
% compute relevant quantitites
n = zeros(s,1); m = n;
for t=1:s
pidx = idx == u(t);
n(t) = sum(pidx);
d = circ_dist(alpha(pidx),med);
m(t) = sum(d<0);
end
if any(n<10)
warning('Test not applicable. Sample size in at least one group to small.') %#ok<WNTAG>
end
M = sum(m);
P = (N^2/(M*(N-M))) * sum(m.^2 ./ n) - N*M/(N-M);
pval = 1 - chi2cdf(P,s-1);
if pval < 0.05
med = NaN;
end
function [alpha, idx] = processInput(varargin)
if nargin==2 && sum(abs(round(varargin{2})-varargin{2}))>1e-5
alpha1 = varargin{1}(:);
alpha2 = varargin{2}(:);
alpha = [alpha1; alpha2];
idx = [ones(size(alpha1)); 2*ones(size(alpha2))];
elseif nargin==2
alpha = varargin{1}(:);
idx = varargin{2}(:);
if ~(size(idx,1)==size(alpha,1))
error('Input dimensions do not match.')
end
else
error('Invalid use of circ_wwtest. Type help circ_wwtest.')
end