forked from circstat/circstat-matlab
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcirc_raotest.m
130 lines (109 loc) · 4.04 KB
/
circ_raotest.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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
function [p U UC] = circ_raotest(alpha)
% [p U UC] = circ_raotest(alpha)
% Calculates Rao's spacing test by comparing distances between points on
% a circle to those expected from a uniform distribution.
%
% H0: Data is distributed uniformly around the circle.
% H1: Data is not uniformly distributed around the circle.
%
% Alternative to the Rayleigh test and the Omnibus test. Less powerful
% than the Rayleigh test when the distribution is unimodal on a global
% scale but uniform locally.
%
% Due to the complexity of the distributioin of the test statistic, we
% resort to the tables published by
% Russell, Gerald S. and Levitin, Daniel J.(1995)
% 'An expanded table of probability values for rao's spacing test'
% Communications in Statistics - Simulation and Computation
% Therefore the reported p-value is the smallest alpha level at which the
% test would still be significant. If the test is not significant at the
% alpha=0.1 level, we return the critical value for alpha = 0.05 and p =
% 0.5.
%
% Input:
% alpha sample of angles
%
% Output:
% p smallest p-value at which test would be significant
% U computed value of the test-statistic u
% UC critical value of the test statistic at sig-level
%
%
% References:
% Batschelet, 1981, Sec 4.6
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
alpha = alpha(:);
% for the purpose of the test, convert to angles
alpha = circ_rad2ang(alpha);
n = length(alpha);
alpha = sort(alpha);
% compute test statistic
U = 0;
lambda = 360/n;
for j = 1:n-1
ti = alpha(j+1) - alpha(j);
U = U + abs(ti - lambda);
end
tn = (360 - alpha(n) + alpha(1));
U = U + abs(tn-lambda);
U = (1/2)*U;
% get critical value from table
[p UC] = getVal(n,U);
function [p UC] = getVal(N, U)
% Table II from Russel and Levitin, 1995
alpha = [0.001, .01, .05, .10];
table = [ 4 247.32, 221.14, 186.45, 168.02;
5 245.19, 211.93, 183.44, 168.66;
6 236.81, 206.79, 180.65, 166.30;
7 229.46, 202.55, 177.83, 165.05;
8 224.41, 198.46, 175.68, 163.56;
9 219.52, 195.27, 173.68, 162.36;
10 215.44, 192.37, 171.98, 161.23;
11 211.87, 189.88, 170.45, 160.24;
12 208.69, 187.66, 169.09, 159.33;
13 205.87, 185.68, 167.87, 158.50;
14 203.33, 183.90, 166.76, 157.75;
15 201.04, 182.28, 165.75, 157.06;
16 198.96, 180.81, 164.83, 156.43;
17 197.05, 179.46, 163.98, 155.84;
18 195.29, 178.22, 163.20, 155.29;
19 193.67, 177.08, 162.47, 154.78;
20 192.17, 176.01, 161.79, 154.31;
21 190.78, 175.02, 161.16, 153.86;
22 189.47, 174.10, 160.56, 153.44;
23 188.25, 173.23, 160.01, 153.05;
24 187.11, 172.41, 159.48, 152.68;
25 186.03, 171.64, 158.99, 152.32;
26 185.01, 170.92, 158.52, 151.99;
27 184.05, 170.23, 158.07, 151.67;
28 183.14, 169.58, 157.65, 151.37;
29 182.28, 168.96, 157.25, 151.08;
30 181.45, 168.38, 156.87, 150.80;
35 177.88, 165.81, 155.19, 149.59;
40 174.99, 163.73, 153.82, 148.60;
45 172.58, 162.00, 152.68, 147.76;
50 170.54, 160.53, 151.70, 147.05;
75 163.60, 155.49, 148.34, 144.56;
100 159.45, 152.46, 146.29, 143.03;
150 154.51, 148.84, 143.83, 141.18;
200 151.56, 146.67, 142.35, 140.06;
300 148.06, 144.09, 140.57, 138.71;
400 145.96, 142.54, 139.50, 137.89;
500 144.54, 141.48, 138.77, 137.33;
600 143.48, 140.70, 138.23, 136.91;
700 142.66, 140.09, 137.80, 136.59;
800 142.00, 139.60, 137.46, 136.33;
900 141.45, 139.19, 137.18, 136.11;
1000 140.99, 138.84, 136.94, 135.92 ];
ridx = find(table(:,1)>=N,1);
cidx = find(table(ridx,2:end)<U,1);
if ~isempty(cidx)
UC = table(ridx,cidx+1);
p = alpha(cidx);
else
UC = table(ridx,end-1);
p = .5;
end