Circular Statistics Toolbox (Directional Statistics)

Compute descriptive and inferential statistics for circular or directional data.

Description

CircStat for Matlab
=======================

Toolbox for circular statistics with Matlab.

Authors: Philipp Berens
Email: [e-mail]
Homepage: http://philippberens.wordpress.com/code/circstats/

Contributors:
Marc Velasco, Tal Krasovsky

Reference:
P. Berens, CircStat: A Matlab Toolbox for Circular Statistics, Journal of Statistical Software, Volume 31, Issue 10, 2009
http://www.jstatsoft.org/v31/i10

Please cite this paper when the provided code is used. See licensing terms for details.

Contents:
circ_r Resultant vector length
circ_mean Mean direction of a sample of circular data
circ_axial Mean direction for axial data
circ_median Median direction of a sample of circular data
circ_std Dispersion around the mean direction (std, mardia)
circ_var Circular variance
circ_skewness Circular skewness
circ_kurtosis Circular kurtosis
circ_moment Circular p-th moment
circ_dist Distances around a circle
circ_dist2 Pairwise distances around a circle
circ_confmean Confidence intervals for mean direction
circ_stats Summary statistics

circ_rtest Rayleigh's test for nonuniformity
circ_otest Hodges-Ajne test (omnibus test) for nonuniformity
circ_raotest Rao's spacing test for nonuniformity
circ_vtest V-Test for nonuniformity with known mean direction
circ_medtest Test for median angle
circ_mtest One-sample test for specified mean direction
circ_wwtest Multi-sample test for equal means, one-factor ANOVA
circ_hktest Two-factor ANOVA
circ_ktest Test for equal concentration parameter
circ_symtest Test for symmetry around median angle
circ_kuipertest Test whether two distributions are identical (like KS test)

circ_corrcc Circular-circular correlation coefficient
circ_corrcl Circular-linear correlation coefficient

circ_kappa Compute concentration parameter of a vm distribution

circ_plot Visualization for circular data
circ_clust Simple clustering for circular data
circ_samplecdf Evaluate CDF of a sample of angles

rad2ang Convert radian to angular values
ang2rad Convert angular to radian values

All functions take arguments in radians (expect for ang2rad). For a detailed description of arguments and outputs consult the help text in the files.

Since 2010, most functions for descriptive statistics can be used in Matlab style matrix computations. As a last argument, add the dimension along which you want to average. This changes the behavior slightly from previous relaeses, in that input is not reshaped anymore into vector format. Per default, all computations are performed columnwise (along dimension 1). If you prefer to use the old functions, for now they are contained in the subdirectory 'old'.

References:
- E. Batschelet, Circular Statistics in Biology, Academic Press, 1981
- N.I. Fisher, Statistical analysis of circular data, Cambridge University Press, 1996
- S.R. Jammalamadaka et al., Topics in circular statistics, World Scientific, 2001
- J.H. Zar, Biostatistical Analysis, Prentice Hall, 1999

If you have suggestions, bugs or feature requests or want to contribute code, please email me.

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Circular Cross Correlation
This submission has inspired the following:
Kernel smoothing density estimate for circular data

• Required Products: Statistics Toolbox
• MATLAB release: MATLAB 7.10 (2010a)

[sanjay sane]

thanks.

[Sooho Park]

Exactly what I needed

[Lyle Muller]

Would include skewness and kurtosis!

[Adrian Bartlett]

I believe I have found an error in one of your functions:
circ_rtest.m
Line 51: z = R^2 / n;
should be
Line 51: z = R^2 * n;

[Adrian B]

I take that back -- guilty of confusing 'r' and 'R'

[Lorenzo Guerrasio]

I wish I saw it before :)
Very nicely done.

[Lorenzo Guerrasio]

PS:
I think there is an error in the circ_dist function

I think this

r = angle(repmat(exp(1i*x(:)'),length(y),1) ...
       ./ repmat(exp(1i*y(:)),1,length(x)));

Let me know if it's correct

[Richard Heitz]

this is a great toolbox

[Tal Krasovsky]

Excellent toolbox, helped me a lot. Greatly appreciated!

[Richard Courtemanche]

This is a great submission, filling an obvious gap in the statistical world out there. Easy to use, well done, and the author provides great feedback.

Great!

[Richard Heitz]

Sorry, but I think circ_var returns s = (1-r) when it should be s = 2*(1-r).

for this reason, circ_std and circ_var will not agree

[Philipp Berens]

Both definitions are around... I will optionally add computing both with the next update.

[Flurin Honegger]

Good work! So far it helped me a lot. But there are some
errors arround!

1.) function circ_hktest

Line 55 found -> qm = zeros(p,1); qr = qm; qn = pm;
   corrected? -> qm = zeros(p,1); qr = qm; qn = qm;

Line 94 found -> eff_2 = sum(qr.^2 ./ sum(cn,2)) - tr.^2/n;
   corrected? -> eff_2 = sum(qr.^2 ./ sum(cn,1)') - tr.^2/n;

Line 107 found -> beta = 1/(1-1/(5*kk)-1/(10*(kk^2)));
       comment -> beta overloads the beta function (help beta)
                  An other name like betaF should be used

Line 144 found -> F1 = beta * ms_1 / ms_r;
       comment -> if inter is set to 0/false beta is not defined!
       
       
2.) function circ_std

The documentation in the paper

  In CircStat , theangular deviation is computed as
  >>s=circ_std(alpha);
  and the circular standard deviation as
  >>s0=circ_std(alpha,[],[],'mardia');
  
does not fit with neither - paper and Matlab Central - implementation.

[Philipp Berens]

Dear Florin,

thanks for the error report.

With regards to 1: Fixed all bugs. I tested the output on the example in Harrison & Kanji.

With regards to 2: This is unfortunate. The current (and more recent) toolbox version returns both, angular deviation and circular standard deviation as first and second return argument.

Appreciate your feedback,

Philipp

[Shiquan Wang]

This is great work! Thanks.

[Shiquan Wang]

bug report:
function stats = circ_stats(alpha, w, d)
line50:stats.std_mardia = circ_std(alpha,w,d,'mardia');
the function circ_std(alpha, w, d) doesn't accept parameter 'mardia'.

[omzaz]

Encountering same issue as Shiquan

[Philipp Berens]

Issue is fixed in the upload of 11/5/09.

[chairmanK]

Functions do not gracefully handle N-dimensional arrays for powerful MATLAB-style computations; instead, inputs are coerced to be column vectors. There are also numerous errors. One example, in circ_moment.m:
cbar = sum(cos(p*alpha'*w))/n;
(p*alpha'*w) is a SCALAR dot product, so clearly this is not a weighted sum of cosines as it ought to be.
There are many other bugs like this. Please fix!

[Cesare]

something's wrong in the new circ_vmpdf....can't replicate results from previous release

[Tomo]

I think, in circ_skewness, the eponent for the denominator of the last equation should be (3/2), rather than (2/3), according to formula (2.29) of Fisher 1993.

[Philipp Berens]

I fixed the bugs reported in December in the first upload of 2010.

I also added the functionality asked for by chairmanK. The functions for descriptive statistics now handle N dimensional arrays and the computations can be performed 'Matlab-style'. As a backup, the new release comes with a folder 'old', which contains the functions that are thus replaced as backup. If you experience problems or issues with the new functions or would like to see additional functions converted let me know. Unused arguments in between can be left empty.

circ_vmpdf seems to work fine with me and produces data with the correct moments. Please be more specific.

[Cesare]

also with the 2010b.
Thanks,
Cesare

[Cesare]

Ops sorry I must have meesed up with the posts....
I'll repost my doubt properly:

Try to run

points = -pi:((4*pi)/(2*Nbin)):3*pi;
mu = -2.838;
kappa = 0.5125;

p = circ_vmpdf(points(1:end-1).', mu, kappa);

It seems to me that results from version circStat2009 (which I think were correct) differ from those of CircStat2009d and CircStat2010b. Maybe I'm doing something wrong. Please let me know.
Thanks a lot,
Cesare

[Philipp Berens]

There has been a slight (and unfortunately undocumented) change in semantics from 2009 to the later versions.

vmpdf computes the density, i.e. it evaluates the probability density function of the von mises distribution at the designated points. The earlier version computed the approximate probability in a small bin with width (alpha(2)-alpha(1) ), as is needed if you want to plot histograms. As you will see, you can easily recover the old behavior by

p = circ_vmpdf(points(1:end-1).', mu, kappa);
p = p * diff(points(1:2));

to obtain approximate probabilities.

[Aisling]

This toolbox is great!! I'm getting an error with circ_clust though. Is there a bug? I haven't been able to solve it myself. This is the error message i get.

??? In an assignment A(I) = B, the number of elements in B and
 I must be the same.

Error in ==> circ_clust at 53
            mu(j) = circ_mean(alpha(cid==j));

[Philipp Berens]

I fixed the bug in circ_clust.

[Kourosh]

Excellent toolbox!!! Thanks Philipp!

[Carlos]

Nice work

[Christopher]

I found circ_median very slow. There might be a more efficient algorithm to get the median, but at least note that around line 42 m1 and m2 are each determined by calculating the same circ_dist2(beta,beta). circ_dist2 takes a LONG time. I would use an intermediate variable, and calculate circ_dist2 only once (to almost halve the time for the function to run, down from 2 seconds to 1 second for 2000 data points).

[Christopher]

I am seeing an anomaly in circ_plot.m that does not make sense to me, can someone please explain or concur that it is a bug.

According to the help hist can plot either count or normalized bins.
In line 110 of circ_plot rose is called to calculate the bins
110] [t,r] = rose(alpha,x);
and in line 112 the normalized bins are plotted:
112] polar(t,r/sum(r),formats)
Now the vectors t and r returned by rose are such that they can be used to plot the bins directly and have the layout [0 n1 n1 0 0 n2 n2 0 0 ...] i.o.w. each element of bin count appears twice, and sum(r) is equal to 2.*length(alpha). So to truly normalize the r, we should divide by half of sum(r) and each bin should be twice as tall. Check by comparison to hist, which returns the same kind of information, but for the first and last bins.
I propose the following replacements for lines 112 and 113
112] polar(t,2.*r./sum(r),formats)
113] mr = max(2.*r./sum(r));

[Christopher]

Another issue with circ_plot.m:
A typo in line 121, should read
      s = varargin{3};
(instead of vargin{1})

[omzaz]

 

[Philipp Berens]

Hi Christopher, thanks for your feedback. I will update circ_plot with the next upload.

[Bart Geurten]

Needed a circular statistic means and got what I needed... thanks!

[Quang Nguyen Dinh]

I'm sorry for this silly question. When I use circ_vmpdf to calculate for the case kappa is big and data has close value with the mean, it return the value bigger than 1. For example:

circ_vmpdf(pi,pi,35)
ans =
    2.3516

Did I mis-understand anything here?
I saw in the graphs of two functions based on kappa:

f1 = 1/(2*pi*besseli(0,kappa))
f2 = exp(kappa)

With high value of kappa, the increasing rate of second one is much higher than the first one so that it's not strange if the above case happened.

[Christopher]

Sorry, the error is in circ_kuipertest.m, not circ_kuiper.m as previously written!

[Christopher]

Um, what happened to my comment? What I wrote was that an typo error appears to have been introduced in circ_kuipertest.m in the advance to version 2011f. In 2010e, line 48 of the file reads:
[phis2 cdf2 phiplot2 cdfplot2] = circ_samplecdf(alpha2, res);

and in version 2011f, that line reads:
[~, cdf2 phiplot2 cdfplot2] = circ_samplecdf(alpha2, res);

and matlab complains of incorrect statement or expression.

[Heida Maria Sigurdardottir]

Thanks people, this toolbox is really helpful and easy to use.

I have one stats question -- forgive me that it is not a direct question about this toolbox but perhaps someone could help nonetheless.

I have repeated measures of circular data for multiple participants, lets say 15 participants where each participant contributes four angles. I have reason to believe that the angular distributions are going to be multipolar and not von Mises distributed. This would in essence require some kind of non-parametric repeated measures test which I am not sure has been developed for circular data. Is there a way to test for circular uniformity in this data set by using the procedures in the circstat toolbox, perhaps using some kind of p-value correction?

[Vlad Atanasiu]

Good toolbox! I added a function for kernel smoothing density estimate for circular data here: http://www.mathworks.com/matlabcentral/fileexchange/32614-kernel-smoothing-density-estimate-for-circular-data .

[omzaz]

Very useful toolbox. Option to ignore NaNs in the calculations would make it even better.

[omzaz]

Can any of the multi-sample tests in this toolbox be used with repeated measures data or do they all assume independent samples?

[Philipp Berens]

Thanks for the comments.

@Christopher: The ~ has been introduced as a placeholder in the latest MATLAB versions for output arguments that are not needed. I will go back to some dummy variable with the next upload.

@Heida: I don't see an easy way of doing what you suggest with the functions implemented.

@Omzaz: The multi-sample tests assume independent samples. I don't know about repeated-measures ANOVA etc. for circular data. If you find anything let me know.

The option to ignore NaNs... I think this is a tricky thing, because you always make a specific choice how NaNs are treated and each user might have different preferences. I will think about it though.

[Fuh-Cherng]

I am new in circular statistics, so don't laugh at me... But I do have a question about the circ_median() function.

Say I have a data set that contains six angles [0.1 0.2 0.3 0.4 0.5 0.6]. when I feed these data into circ_median(), the function returns a median = 0.4

I thought that, when a data set contains an even number of observations, the median would be calculated as the average of the middle two numbers (i.e., (0.3+0.4)/2 = 0.35).

My code is listed below.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = [0.1 0.2 0.3 0.4 0.5 0.6]';
med = circ_median(alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Can anyone help me with this?

Sincerely,
Fuh

[Christopher]

@Fuh: indeed it should and when I step carefully through the function, sometimes the result comes out correct and sometimes it doesn't, somewhat dependent on the numbers in alpha. To fix the problem go to lines 45 and 46 of circ_median (ver 2011f). You see two inequalities, dd>=0 and dd<0. The two inequalities should be identical for consistency and the correct result. Edit line 46 to read:
    m2 = sum(dd<=0,1);
Now the function seems to behave as expected.

[Fuh-Cherng]

@Christopher: Thank you so much for your kindness and help. I really appreciate it.

[Balazs Barkoczi]

Thanks for this excellent toolbox!
I have only some problems with the example files, that I downloaded from http://www.jstatsoft.org/v31/i10
example1:
??? Undefined function or method 'parseVarArgs', therefore the figure 2 isn't complete, and it hasn't axis labels.

example2:
??? Error using ==> mtimes
Inner matrix dimensions must agree.

Error in ==> example2 at 42
  zm = r*exp(i*phi);
Perhaps a dot is absent, but after this modification zm = r.*exp(i*phi); the same error occurs:
??? Undefined function or method 'parseVarArgs'
Can somebody help me to fix this problems?
Thank you very much!

[Luke]

Great tool.
I do have to say that circ_mtest is a bit weird.
The input is [pval, z] but output is set as [h,mu,ul,ll]