Skip to content

Commit b8ea0b0

Browse files
dervinismdervinism
authored
Add files via upload
1 parent a2ee73c commit b8ea0b0

3 files changed

Lines changed: 201 additions & 0 deletions

File tree

‎WatsonU2_goodnessOfFit.m‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
function [U2, pval] = WatsonU2_goodnessOfFit(F, sumTerms)
2+
% [U2, pval] = WatsonU2_goodnessOfFit(F, sumTerms)
3+
%
4+
% Function computes Watson U2 goodness-of-fir statistic.
5+
% Input: F - fitted cummulative distribution fucntion of your fit.
6+
% sumTerms - number of exponential sum terms in the estimation of
7+
% p-value. The default value is 30 as a larger number of
8+
% terms provides negligible improvement in convergence.
9+
% Output: U2 - Watson U2 goodness-of-fit statistic.
10+
% pval - associated p-value. Larger U2 values give smaller
11+
% p-values. Therefore, the bigger is the U2 statistic, the
12+
% poorer the fit. If your aim is to test the similarity of
13+
% the fitted and empirical distributions, you should,
14+
% therefore, invert your p-value:
15+
% new p-value = 1 - original p-value.
16+
17+
% by Martynas Dervinis (martynas.dervinis@gmail.com)
18+
19+
if nargin < 2
20+
sumTerms = 30;
21+
end
22+
23+
n = numel(F);
24+
u_bar = mean(F);
25+
26+
% Rank the sample
27+
[~, sortOrder] = sort(F);
28+
ranks = zeros(1,n);
29+
for iF = 1:n
30+
ranks(sortOrder(iF)) = iF;
31+
end
32+
33+
% Calculate U2 goodness-of-fit statistic
34+
U2 = 0;
35+
for iF = 1:n
36+
U2 = U2 + (F(iF) - (u_bar - 0.5) - ((2*ranks(iF) - 1)/(2*n)))^2;
37+
end
38+
U2 = U2 + (1/(12*n));
39+
40+
% Estimate the p-value
41+
pval = 0;
42+
for iSum = 1:sumTerms
43+
pval = pval + ((-1)^(iSum-1))*2*exp(-2*(iSum^2)*(pi^2)*U2);
44+
end

‎circ_ksdensity.m‎

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
function [vfEstimate] = circ_ksdensity(vfObservations, vfPDFSamples, vfDomain, fSigma, vfWeights)
2+
3+
% circ_ksdensity FUNCTION - Compute a kernel density estimate over a periodic domain
4+
%
5+
% Usage: [vfEstimate] = circ_ksdensity(vfObservations, vfPDFSamples,
6+
% <vfDomain, fSigma, vfWeights>)
7+
%
8+
% This function calculates a kernel density estimate of an (optionally
9+
% weighted) data sample, over a periodic domain.
10+
%
11+
% 'vfObservations' is a set of observations made over a periodic domain,
12+
% optionally defined by 'vfDomain': [fMin fMax]. The default domain is
13+
% [0..2*pi]. 'vfPDFSamples' defines the sample points over which to perform
14+
% the kernel density estimate, over the same domain as 'vfObservations'.
15+
%
16+
% Weighted estimations can be performed by providing the optional argument
17+
% 'vfWeights', where each element in 'vfWeights' corresponds to the
18+
% matching element in 'vfObservations'.
19+
%
20+
% The kernel density estimate will be performed using a wrapped Gaussian
21+
% kernel, with a width estimated as
22+
% (4/3)^0.2 * circ_std(vfObservations, vfWeights) * (length(vfObservations^-0.2)
23+
%
24+
% The optional argument 'fSigma' can be provided to set the width of the
25+
% kernel.
26+
%
27+
% 'vfEstimate' will be a vector with a (weighted) estimate of the
28+
% underlying distribution, with an entry for each element of
29+
% 'vfPDFSamples'. If no weighting is supplied, the estimate will be scaled
30+
% such that it forms a PDF estimate over the supplied sample domain, taking
31+
% into account sample bin widths. If a weight vector is supplied then the
32+
% estimate will be scaled such that the sum over the domain attempts to
33+
% match the sum of weights, taking into account sample bin widths.
34+
35+
% Author: Dylan Muir <dylan.muir@unibas.ch>
36+
% Created: 23rd October, 2013
37+
38+
% -- Defaults
39+
40+
DEF_vfDomain = [0 2*pi];
41+
42+
43+
% -- Check arguments
44+
45+
if (nargin < 2)
46+
help circ_ksdensity;
47+
error('circ_ksdensity:Usage', ...
48+
'*** circ_ksdensity: Incorrect usage');
49+
end
50+
51+
if (~exist('vfDomain', 'var') || isempty(vfDomain))
52+
vfDomain = DEF_vfDomain;
53+
end
54+
55+
% - Do we need to estimate fSigma?
56+
if (~exist('fSigma', 'var'))
57+
% - Sigma will be estimated
58+
fSigma = [];
59+
end
60+
61+
vfObservations = vfObservations(:);
62+
vnPSFSamplesSize = size(vfPDFSamples);
63+
vfPDFSamples = vfPDFSamples(:);
64+
65+
% - If weights are not provided, weight each observation equally
66+
if (~exist('vfWeights', 'var'))
67+
vfWeights = ones(size(vfObservations)) ./ numel(vfObservations);
68+
69+
% - Check the number of observations matches the number of weights
70+
elseif (numel(vfObservations) ~= numel(vfWeights))
71+
error('circ_ksdensity:Usage', ...
72+
'*** circ_ksdensity: The number of observations must be equal to the number of weights.');
73+
end
74+
75+
76+
% -- Map everything to [0..2pi] and wrap over domain
77+
78+
vfObservations = (vfObservations - vfDomain(1)) ./ diff(vfDomain) .* 2*pi;
79+
vfObservations = mod(vfObservations, 2*pi);
80+
81+
vfPDFSamples = (vfPDFSamples - vfDomain(1)) ./ diff(vfDomain) .* 2*pi;
82+
vfPDFSamples = mod(vfPDFSamples, 2*pi);
83+
84+
85+
% - Estimate sigma, if necessary
86+
if (isempty(fSigma))
87+
fSigma = (4/3)^0.2 * circ_std(vfObservations, vfWeights) * (numel(vfObservations)^-0.2);
88+
end
89+
90+
91+
% -- Pad observations above and below domain
92+
93+
vfObservations = [vfObservations;
94+
vfObservations - 2*pi;
95+
vfObservations + 2*pi];
96+
vfWeights = repmat(vfWeights, 3, 1) ./ 3;
97+
98+
99+
% -- Perform kernel density estimate
100+
101+
vfEstimate = ksdensity(vfObservations, vfPDFSamples, 'weights', vfWeights', 'width', fSigma);
102+
103+
% - Reshape return to match shape of 'vfPDFSamples'
104+
vfEstimate = reshape(vfEstimate, vnPSFSamplesSize);
105+
106+
% - Correct scaling of histogram estimate
107+
vfEstimate = vfEstimate .* 3 .* sum(vfWeights);
108+
109+
% --- END of circ_ksdensity FUNCTION ---
110+
111+
112+
% circ_std FUNCTION Estimate the weighted circular standard deviation of a dataset
113+
function [s s0] = circ_std(alpha, w)
114+
115+
% compute mean resultant vector length
116+
r = circ_r(alpha,w);
117+
118+
s = sqrt(2*(1-r)); % 26.20
119+
s0 = sqrt(-2*log(r)); % 26.21
120+
121+
122+
% circ_r FUNCTION Compute the weighted resultant of a dataset
123+
function r = circ_r(alpha, w)
124+
125+
% compute weighted sum of cos and sin of angles
126+
r = sum(w.*exp(1i*alpha),1);
127+
128+
% obtain length
129+
r = abs(r)./sum(w,1);
130+
131+
132+
% --- END of circ_ksdensity.m ---
133+
134+

‎recentrePhase.m‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
function phaseData = recentrePhase(phaseData, phaseCentre)
2+
% phaseData = recentrePhase(phaseData, phaseCentre)
3+
%
4+
% Function converts phase values to values within a chosen phase range
5+
% centred around input variable phaseCentre: [phaseCentre-pi
6+
% phaseCentre+pi].
7+
% Input: phaseData - a vector or a matrix of phase values.
8+
% phaseCentre - the centre of the new phase range.
9+
% Output: phaseData - recentred phase data.
10+
11+
upperLimit = phaseCentre + pi;
12+
lowerLimit = phaseCentre - pi;
13+
for i = 1:numel(phaseData)
14+
if phaseData(i) > upperLimit
15+
while phaseData(i) > upperLimit
16+
phaseData(i) = phaseData(i) - 2*pi;
17+
end
18+
elseif phaseData(i) < lowerLimit
19+
while phaseData(i) < lowerLimit
20+
phaseData(i) = phaseData(i) + 2*pi;
21+
end
22+
end
23+
end

0 commit comments

Comments
 (0)