forked from circstat/circstat-matlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
circ_samplecdf.m
74 lines (55 loc) · 1.69 KB
/
circ_samplecdf.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
function [phis, cdf, phiplot, cdfplot] = circ_samplecdf(thetas, resolution)
% [phis, cdf, phiplot, cdfplot] = circ_samplecdf(thetas, resolution)
%
% Helper function for circ_kuipertest.
% Evaluates CDF of sample in thetas.
%
% Input:
% thetas sample (in radians)
% resolution resolution at which the cdf is evaluated
%
% Output:
% phis angles at which CDF is evaluated
% cdf CDF values at these angles
% phiplot as phi, for plotting
% cdfplot as cdf, for plotting
%
%
% Circular Statistics Toolbox for Matlab
% By Marc J. Velasco, 2009
if nargin < 2
resolution = 100;
end
phis = 0;
cdf = zeros(1, length(phis));
phis = linspace(0,2*pi,resolution+1);
phis = phis(1:end-1);
% ensure all points in thetas are on interval [0, 2pi)
x = thetas(thetas<0);
thetas(thetas<0) = (2*pi-abs(x));
% compute cdf
thetas = sort(thetas);
dprob = 1/length(thetas); %incremental change in probability
cumprob = 0; %cumultive probability so far
% for a little bit, we'll add on 2pi to the end of phis
phis = [phis 2*pi];
for j=1:resolution
minang = phis(j);
maxang = phis(j+1);
currcount = sum(thetas >= minang & thetas < maxang);
cdf(j) = cumprob + dprob*currcount;
cumprob = cdf(j);
end
phis = phis(1:end-1);
% for each point in x, duplicate it with the preceding value in y
phis2 = phis;
cdf2 = [0 cdf(1:end-1)];
cdfplottable = [];
phisplottable = [];
for j=1:length(phis);
phisplottable = [phisplottable phis(j) phis2(j)]; %#ok<AGROW>
cdfplottable = [cdfplottable cdf2(j) cdf(j)]; %#ok<AGROW>
end
phiplot = [phisplottable 2*pi];
cdfplot = [cdfplottable 1];