forked from circstat/circstat-matlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
circ_otest.m
81 lines (61 loc) · 1.85 KB
/
circ_otest.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
function [pval m] = circ_otest(alpha, sz, w)
%
% [pval, m] = circ_otest(alpha,sz,w)
% Computes Omnibus or Hodges-Ajne test for non-uniformity of circular data.
% H0: the population is uniformly distributed around the circle
% HA: the population is not distributed uniformly around the circle
%
% Alternative to the Rayleigh and Rao's test. Works well for unimodal,
% bimodal or multimodal data. If requirements of the Rayleigh test are
% met, the latter is more powerful.
%
% Input:
% alpha sample of angles in radians
% [sz step size for evaluating distribution, default 1 degree
% [w number of incidences in case of binned angle data]
% Output:
% pval p-value
% m minimum number of samples falling in one half of the circle
%
% PHB 3/16/2009
%
% References:
% Biostatistical Analysis, J. H. Zar
% A bivariate sign test, J. L. Hodges et al., 1955
% A simple test for uniformity of a circular distribution, B. Ajne, 1968
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% [email protected] - www.kyb.mpg.de/~berens/circStat.html
if size(alpha,2) > size(alpha,1)
alpha = alpha';
end
if nargin < 2 || isempty(sz)
sz = circ_ang2rad(1);
end
if nargin < 3
w = ones(size(alpha));
else
if length(alpha)~=length(w)
error('Input length does not match.')
end
w =w(:);
end
alpha = mod(alpha,2*pi);
n = sum(w);
dg = 0:sz:pi;
m1 = zeros(size(dg));
m2 = zeros(size(dg));
for i=1:length(dg)
m1(i) = sum((alpha > dg(i) & alpha < pi + dg(i)).*w);
m2(i) = n - m1(i);
end
m = min(min([m1;m2]));
if n > 50
% approximation by Ajne (1968)
A = pi*sqrt(n) / 2 / (n-2*m);
pval = sqrt(2*pi) / A * exp(-pi^2/8/A^2);
else
% exact formula by Hodges (1955)
pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);
end