Code covered by the BSD License

# Circular Statistics Toolbox (Directional Statistics)

### Philipp Berens (view profile)

08 Apr 2006 (Updated )

Compute descriptive and inferential statistics for circular or directional data.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

circ_hktest(alpha, idp, idq, inter, fn)
```function [pval table] = circ_hktest(alpha, idp, idq, inter, fn)

%
% [pval, stats] = circ_hktest(alpha, idp, idq, inter, fn)
%   Parametric two-way ANOVA for circular data with interations.
%
%   Input:
%     idp     indicates the level of factor 1 (1:p)
%     idq     indicates the level of factor 2 (1:q)
%     inter   0 or 1 - whether to include effect of interaction or not
%     fn      cell array containing strings with the names of the factors
%
%
%   Output:
%     pval    vector of pvalues testing column, row and interaction effects
%     table   cell array containg the anova table
%
%   The test assumes underlying von-Mises distributrions.
%   All groups are assumed to have a common concentration parameter k,
%   between 0 and 2.
%
% Update 2012
% PHB 7/19/2009 with code by Tal Krasovsky, Mc Gill University
%
% References:
%   Harrison, D. and Kanji, G. K. (1988). The development of analysis of variance for
%   circular data. Journal of applied statistics, 15(2), 197-223.
%
% Circular Statistics Toolbox for Matlab

% process inputs
alpha = alpha(:); idp = idp(:); idq = idq(:);

if nargin < 4
inter = true;
end

if nargin < 5
fn = {'A','B'};
end

% number of groups for every factor
pu = unique(idp);
p = length(pu);
qu = unique(idq);
q = length(qu);

% number of samples
n = length(alpha);

% compute important sums for the test statistics
cn = zeros(p,q); cr = cn;
pm = zeros(p,1); pr = pm; pn = pm;
qm = zeros(q,1); qr = qm; qn = qm;
for pp = 1:p
p_id = idp == pu(pp); % indices of factor1 = pp
for qq = 1:q
q_id = idq == qu(qq); % indices of factor2 = qq
idx = p_id & q_id;
cn(pp,qq) = sum(idx);     % number of items in cell
cr(pp,qq) = cn(pp,qq) * circ_r(alpha(idx)); % R of cell
end
% R and mean angle for factor 1
pr(pp) = sum(p_id) * circ_r(alpha(p_id));
pm(pp) = circ_mean(alpha(p_id));
pn(pp) = sum(p_id);
end

% R and mean angle for factor 2
for qq = 1:q
q_id = idq == qu(qq);
qr(qq) = sum(q_id) * circ_r(alpha(q_id));
qm(qq) = circ_mean(alpha(q_id));
qn(qq) = sum(q_id);
end

% R and mean angle for whole sample (total)
tr = n * circ_r(alpha);

% estimate kappa
kk = circ_kappa(tr/n);

% different formulas for different width of the distribution
if kk > 2
% large kappa

% effect of factor 1
eff_1 = sum(pr.^2 ./ sum(cn,2)) - tr.^2/n;
df_1 = p-1;
ms_1 = eff_1 / df_1;

% effect of factor 2
eff_2 = sum(qr.^2 ./ sum(cn,1)') - tr.^2/n;
df_2 = q-1;
ms_2 = eff_2 / df_2;

% total effect
eff_t = n - tr.^2/n;
df_t = n-1;

m = mean(cn(:));

if inter

% correction factor for improved F statistic
beta = 1/(1-1/(5*kk)-1/(10*(kk^2)));

% residual effects
eff_r = n - sum(sum(cr.^2./cn));
df_r = p*q*(m-1);
ms_r = eff_r / df_r;

% interaction effects
eff_i = sum(sum(cr.^2./cn)) - sum(qr.^2./qn) ...
- sum(pr.^2./pn) + tr.^2/n;
df_i = (p-1)*(q-1);
ms_i = eff_i/df_i;

% interaction test statistic
FI = ms_i / ms_r;
pI = 1-fcdf(FI,df_i,df_r);

else

% residual effect
eff_r = n - sum(qr.^2./qn)- sum(pr.^2./pn) + tr.^2/n;
df_r = (p-1)*(q-1);
ms_r = eff_r / df_r;

% interaction effects
eff_i = [];
df_i = [];
ms_i =[];

% interaction test statistic
FI = [];
pI = NaN;
beta = 1;
end

% compute all test statistics as
%  F = beta * MS(A) / MS(R);

F1 = beta * ms_1 / ms_r;
p1 = 1 - fcdf(F1,df_1,df_r);

F2 = beta * ms_2 / ms_r;
p2 = 1 - fcdf(F2,df_2,df_r);

else
% small kappa

% correction factor
rr = besseli(1,kk) / besseli(0,kk);
f = 2/(1-rr^2);

chi1 = f * (sum(pr.^2./pn)- tr.^2/n);
df_1 = 2*(p-1);
p1 = 1 - chi2cdf(chi1, df_1);

chi2 = f * (sum(qr.^2./qn)- tr.^2/n);
df_2 = 2*(q-1);
p2 = 1 - chi2cdf(chi2, df_2);

chiI = f * (sum(sum(cr.^2 ./ cn)) - sum(pr.^2./pn) ...
- sum(qr.^2./qn)+ tr.^2/n);
df_i = (p-1) * (q-1);
% pI = 1 - chi2pdf(chiI, df_i); % bugfix 2012
pI = 1 - chi2cdf(chiI, df_i);

end

na = nargout;
if na < 2
printTable;
end

prepareOutput;

function printTable

if kk>2

fprintf('\nANALYSIS OF VARIANCE TABLE (HIGH KAPPA MODE)\n\n');

fprintf('%s\t\t\t\t%s\t%s\t\t%s\t\t%s\t\t\t%s\n', ' ' ,'d.f.', 'SS', 'MS', 'F', 'P-Value');
fprintf('--------------------------------------------------------------------\n');
fprintf('%s\t\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', fn{1}, df_1 , eff_1, ms_1, F1, p1);
fprintf('%s\t\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', fn{2}, df_2 , eff_2, ms_2, F2, p2);
if (inter)
fprintf('%s\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', 'Interaction', df_i , eff_i, ms_i, FI, pI);
end
fprintf('%s\t\t%u\t\t%.2f\t%.2f\n', 'Residual ', df_r, eff_r, ms_r);
fprintf('--------------------------------------------------------------------\n');
fprintf('%s\t\t%u\t\t%.2f', 'Total   ',df_t,eff_t);
fprintf('\n\n')
else
fprintf('\nANALYSIS OF VARIANCE TABLE (LOW KAPPA MODE)\n\n');

fprintf('%s\t\t\t\t%s\t%s\t\t\t%s\n', ' ' ,'d.f.', 'CHI2', 'P-Value');
fprintf('--------------------------------------------------------------------\n');
fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{1}, df_1 , chi1, p1);
fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{2}, df_2 , chi2, p2);
if (inter)
fprintf('%s\t\t%u\t\t%.2f\t\t\t%.4f\n', 'Interaction', df_i , chiI, pI);
end
fprintf('--------------------------------------------------------------------\n');
fprintf('\n\n')

end

end

function prepareOutput

pval = [p1 p2 pI];

if na > 1
if kk>2
table = {'Source','d.f.','SS','MS','F','P-Value'; ...
fn{1}, df_1 , eff_1, ms_1, F1, p1; ...
fn{2}, df_2 , eff_2, ms_2, F2, p2; ...
'Interaction', df_i , eff_i, ms_i, FI, pI; ...
'Residual', df_r, eff_r, ms_r, [], []; ...
'Total',df_t,eff_t,[],[],[]};
else
table = {'Source','d.f.','CHI2','P-Value'; ...
fn{1}, df_1 , chi1, p1;
fn{2}, df_2 , chi2, p2;
'Interaction', df_i , chiI, pI};
end
end

end

end

```