Editor's Note: This file was selected as MATLAB Central Pick of the Week
CircStat for Matlab
=======================
Toolbox for circular statistics with Matlab.
Authors: Philipp Berens
Email: philipp@bethgelab.org
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 (not the technical report!). 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 pth 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 HodgesAjne test (omnibus test) for nonuniformity
circ_raotest Rao's spacing test for nonuniformity
circ_vtest VTest for nonuniformity with known mean direction
circ_medtest Test for median angle
circ_mtest Onesample test for specified mean direction
circ_wwtest Multisample test for equal means, onefactor ANOVA
circ_hktest Twofactor 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 Circularcircular correlation coefficient
circ_corrcl Circularlinear 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).
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.
1.21.0.0  Update on median, wwtest, hktest, kuiper. Also added examples from paper. 

1.20.0.0  Bugfixes in wwtest, kuiper, median, hktest and added example files from paper. 

1.19.0.0  Updates fixing the bugs reported in the last few months. Touched files:


1.18.0.0  Bugfixes and updates 

1.17.0.0  Bug in circ_clust fixed. 

1.16.0.0  Bugfix in circ_hktest: lines 159 and 163 were switched. 

1.15.0.0  Small bugfixes, mainly in the help sections 

1.14.0.0  Various bugfixes. Added Matlabstyle computations. 

1.13.0.0  Bug fix in circ_stats. 

1.12.0.0  Bug fix in circ_skewness and circ_kurtosis. Thanks to Shiquan Wang. 

1.11.0.0  Bug fix in circ_dist and circ_clust. 

1.10.0.0  Bug fix. 

1.9.0.0  Two new tests 

1.8.0.0  Updated reference for paper 

1.7.0.0  Added reference. Removed some bugs. Added new, more complicated tests (ANOVA like testing). 

1.6.0.0  A number of small bug fixes. 

1.5.0.0  Changed licensing 
Inspired by: Circular Cross Correlation
Inspired: Family of Rayleigh Statistics Toolbox, Alphonse:Rex  Handwriting and Print Expertise, circular corcoeff, Kernel smoothing density estimate for circular data, pierremegevand/watsons_u2, Kernel density estimation for circular functions, Seis_Pick, Local shape outline orientation descriptor, CircHist  circular / polar / angle histogram
Create scripts with code, output, and formatted text in a single executable document.
Mahlega Hassanpour (view profile)
Thanks for sharing. Which function calculates the minor and major eigenvalues?
Alvaro RiveraRei (view profile)
Omer Dauber (view profile)
PAK (view profile)
Hi,
These functions have been extremely helpful, given that I am very new to circular statistics.
One question  I'm trying to calculate the circular mean and 95% CI of a vector of phase values that go from pi to pi. Circ_mean calculates mu fine, but my upper/lower limits return NaNs with the message
Warning: Requirements for confidence levels not met.
> In circ_confmean (line 70)
In circ_mean (line 53)
In Phase_Lock_Recall (line 150)
Would you be able to help with this? Do values need to begin at 0?
Thank you!
PAK
Suddha Sourav (view profile)
This is one of the most wellwritten code I have ever come across. Deeply impressed by the numbering of equations in the code. Thank you very much for the contribution!
Josephine Benjamin (view profile)
Hi i'm new in cthis line of topic but can you help me generate a mixture of vonMises distribution with p=[0.3 0.7], mu=[1.5pi 3.2pi] and kappa=[3.6 2.8] i think circ_vmrnd function only inputs mean directions and concetration parameters. Thanks a lot
Keerthi Nagothu (view profile)
Thank you, I am now able to easily calculate circular median! However I see the below warning:
Warning: Ties detected.
> In circ_median (line 58)
Can you help me understand what it means?
Yu Du (view profile)
These are fantastic codes! The helped me a lot! Thanks!
I have a question about the cluster function circ_clust. I think there is something wrong with it. For example, I tried to use [0 0 0 pi pi pi] to get two clusters. But I did not get cid as [1 1 1 2 2 2]. Instead, I got cid = [1 2 1 1 1 1], which is very weird to me. Could you check the code and add more explanation about it? Thank you in advance.
Vincent BARALE (view profile)
Hi
Thanks for the work, realy usefull and simple to use !!
Daniel de Malmazet (view profile)
Hi,
Thanks for the work! I've noticed a problem with the circ_median function. At line 54, if there are more than 2 elements that are equal to m, you might end up computing an inaccurate mean. It seems better replacing it with: " idx = find(dm==m,n); "
Best
Daniel de Malmazet
Giovanni Esposito (view profile)
hello ,
I have the directions expressed in degrees North but it is possible to use the circ_corr function to calculate the correlation or do you have to go in degrees east?
Thanks in advance
anisse khald (view profile)
hello
can you help me please about formatSubplot function? I can't find her. thank you very much
andrea bertana (view profile)
Sorry, the previous rating was a mistake that I didn't notice until a colleague just pointed out for me!
andrea bertana (view profile)
Abhilash (view profile)
Hey Philipp,
In the circ_vmpdf function, line 45, it currently says this 
C = 1/(2*pi*bessi0(kappa));
I assume it should be besseli(0,kappa) ?
Yavor Kamer (view profile)
@Adrian Bondy you da real MVP
Leslie (view profile)
Very intriguing ... but circ_r and circ_stats (and presumably other functions) return "NaN" if there are any NaN in the input data! This should at the very least be clearly mentioned in the documentation.
Paul (view profile)
Absolutely fantastic, thank you so much for putting all of this in one place.
Alessia Cacace (view profile)
Hello,
I am getting an error with the circ_wwtest that says: "warning: test not applicable. Average resultant vector length less than 0.45 Check assumption line 109"
I am using data points coming from two different rose plots of angular data.
Can anyone help with this?
Thanks!
Nike (view profile)
Philipp Berens (view profile)
@andrea bertana: What did you not like?
@caseyraj: You convert the angle to radians and make a long vector. You need to additional vectors: factor1 and factor 2. In your case, factor1 indicates target and takes values 1, 2 and 3. factor 2 would be condition and takes values 1 and 2.
andrea bertana (view profile)
caseyraj (view profile)
Hello,
I am currently trying to analyze some data consisting of angular data. Therefore I am trying to use the following function from Circular Statistics toolbox:
[p,F] = CircularANOVA(angles,[factor1 factor2],method)
The sample data is provided below as well. Any help on how to use this function would be greatly appreciated
Data:
Target 1 Target 2 Target 3
Condition 1 30.0 40.0 50.0
Condition 1 40.0 50.0 30.0
Condition 2 30.0 34.0 25.0
Condition 2 20.0 34.0 30.0
Philipp Berens (view profile)
Dear all,
please not that feature requests, bug fixes or performance improvements should be made on github at https://github.com/circstat/circstatmatlab.
Thank you!
Tim Kietzmann (view profile)
Perfect toolbox, so helpful. What I am missing is an implementation of the covariance (similar to the "cov" function in matlab). Could this be added by any chance?
Renaud JG (view profile)
I have a question about how to use the function circ_hktest. I want to use it to analyse the difference in the mean of phase relation value between movement of the right vs left leg during locomotion at different time points following a spinal cord injury across different groups. The data comes from several subjects, so when I put all the data in a single column, the degrees of freedom are consistent with the number of groups and timepoints I have in the experimental design put the error match the total number of phas relation value I got, which equals the number of steps taken by all my subjects, which don't make sense. Also, the F is really high because I think the test mix intra and intersubject variability. Does I understand correctly? Thanks for helping me using your great toolbox!
Tristan Chaplin (view profile)
Thanks for the great toolbox.
I have question about the applications to neuroscience. I had low firing cell (3Hz above spontaneous) that looked quite direction selective (CV = 0.3) but it just failed the Rayleigh test (p = 0.07). I notice the magnitude of weights (spike rates) used in the Rayleigh test affect the result  if I multiply the rates by 100 the result is highly significant, so it would seem it might not work well for low firing cells?
Then, I tried using the firing rate of each trial, rather than the mean firing rate for each direction, which your examples seem to use. The result is now high significant.
Do you have any recommendation on handling low firing rates and whether it is more appropriate to use trial rates or mean rates?
Alberto Dallolio (view profile)
Great Job!
I'm trying to compute the standard deviation for angles..
Is there any way to do this with this toolbox?
Thanks a lot!
Alb
Michael Melnychuk (view profile)
Hello, thanks for all your work putting this together. I have a question.
In the script 'circ_wwtest' there are a set of assumptions that are checked. Sample size and average vector length are the parameters. I am wondering what the rationale for this is. Should the Ftest result (say on a very large sample for example) not outweigh arbitrary, hardcoded assumptions?
If you have a reference or a mathematical reason for these assumptions please let me know. :)
Thank you.
Philipp Berens (view profile)
I have added the CircStat toolbox to github and would like you to add all feature requests and bug reports as issues there.
https://github.com/circstat/circstatmatlab
If you have fixes for bugs you can also create pull requests and I will incorporate them.
steve rogers (view profile)
How do I get this toolbox installed in matlab?
Adrian Bondy (view profile)
circ_vmpdf has a serious problem with numerical stability for large kappa (>700 for double precision returns NaN).
This is because the code tries to explicitly evaluate the bessel function in the formula for the von Mises, which is intractable for large kappa. This is unnecessary, since what is ultimately being computed is a ratio between the Bessel function and a large number in the numerator. By calling the bessel function with an additional argument:
besseli(0,kappa,1)
Matlab computes the ratio of the Bessel function and exp(kappa) which is numerically stable.
Then you can simply replace the last two lines with:
p = exp( kappa*(cos(alphathetahat)1)) / (2*pi*besseli(0,kappa,1));
This produces identical results within machine precision and is numerically stable for large kappa.
Shai25 (view profile)
Thanks for the great toolbox.
I have been trying to use the function circ_clust. However, I keep getting the following error message:
In an assignment A(I) = B, the number of elements in B and I must be the same.
Error in circ_clust (line 53)
mu(j) = circ_mean(alpha(cid==j)');
Can you please advise how to overcome this problem?
Chethan Pandarinath (view profile)
C Potgieter (view profile)
Great package.
I found a small inconsistency in the circ_vmrnd program. On line 50, it should read:
alpha = 2*pi*rand(n,1)pi;
Currently the "pi" isn't in the code, which means for large values of kappa, the distribution takes values on the interval (pi,pi), while for kappa close to 0 it takes values between (0,2pi).
Dave (view profile)
This is a great toolbox. However, I have an issue with circ_corrcl, or perhaps I'm misunderstanding its proper use. When I feed it perfectly correlated data, I do not get rho=1:
circ_corrcl(linspace(pi,pi,1000),linspace(0,10,1000))
ans = .7785
I noticed this while trying to see if I could get a sign for the rho value (which is always positive by the definition in the function). Unfortunately, I don't have a copy of the Zar text available.
Owen Brimijoin (view profile)
indispensable
Ch. Lat. (view profile)
you meant "circ_rtest Rayleigh's test for nonuniformity " this is not what you look for ?
Ch. Lat. (view profile)
@ Wasim Malik. I believe the toolbox fieldtrip has it, but their implementation is a bit more bothersome.
Wasim Malik (view profile)
On a quick look, the MooreRayleigh test for uniformity of vector data (B.R. Moore, Biometrika, 1980) does not seem to be available in this toolbox. Philipp, do you have any plans to implement it? Alternatively, does anyone know if a Matlab implementation of that test is available elsewhere? Thanks.
Wasim Malik (view profile)
Diego (view profile)
Hi everybody!
I have a question about circ_plot.m; When I execute this code the angles appear from 0 to 360 degrees.
I only want represent values from 0 to 180. How I can do it? Thanks in advance!
Richard Edmonds (view profile)
Actually, ignore the inverse_cdf function I have provided. It should generate a vlaue for kappa and it needs adjusting for values of thetahat other than zero.
Richard Edmonds (view profile)
Great submission. It would be nice to have cdf and inversion cdf for the vmpdf functions. Here's what I wrote for my needs
function p = circ_vmcdf(alpha, thetahat, kappa)
%integrates the pdf from an angle of pi to an angle alpha
F = @(x)circ_vmpdf(x, thetahat, kappa);
p = quad(F,pi(),alpha);
end
function theta = circ_vminv(p, thetahat, kappa)
%computes the inverse of the abovecirc_vmcdf.
fun =@(alpha)(circ_vmcdf(alpha, thetahat, kappa)p);
theta = fzero(fun,[pi pi]);
end
Sepideh (view profile)
Hi,
Thanks for the great contribution.
Can you please let me know if it is ok to get negative mean or median?
Shall I add 360 to the final angle to make it positive?
Cheers,
Sep
Pete (view profile)
Thanks for the toolkit.
Does anybody have a clue how to do multipleregression with circular data?
bbegginnerr (view profile)
Hi everyone,
Quick question regarding circ_hktest  I quite often get NaNs as an output in the 'Interaction' row. Any idea what am I doing wrong?
Thanks in advance for help.
Philipp Berens (view profile)
Mear  could you be more specific regarding your doubts about the wwtest?
The negative values for circ_mean are a result of the way circ_mean is implemented. If you prefer them to be between 0 and 2pi, just edit the function to provide the data in that format.
I'll update the von Mises function in a future release.
Mear (view profile)
Just a quick question concerning the circ_mean function: I usually get the results in negative values, despite all the input angles being in positive degrees (conv. to rads). It's hardly a big deal to translate this to [0,360] degrees, but it is a bit annoying and seems unneccesary. Is this how it should be? I'm also getting some results for the wwtest which seem very wrong to me (but make sense in light of the negative mean values), and it's making me question the accuracy of this toolbox.
sergio (view profile)
Marc, thanks, you are right.
I generated a von mises distribution with the mu and kappa estimated from my angles, say x, i.e.:
[mu kappa] = circ_vmpar(x)
and then
vonmis = circ_randvm(mu,kappa,length(x))
Then I use the kuiper test to see whether the two distribution x and vonmis differ significantly (the difference can be in any property, such as mean, location and dispersion):
[H,pValue] = circ_kuipertest(x, vonmis)
However I was wondering if it is possible to have more accurate pvalue estimates in the Kuiper test, as already asked by another user before.
sergio (view profile)
Dear Mark,
thanks for your tip, however I'm not really convinced.
Both the circ_ktest and the circ_kuipertest are not described in the pdf:
http://www.jstatsoft.org/v31/i10/paper
Anyway, circ_ktest is a parametric twosample test to determine whether two concentration parameters are different.
The circ_kuipertest is a twosample test which allow to test whether two input samples differ significantly. The difference can be in any property, such as mean location and dispersion. It is a circular analogue of the KolmogorovSmirnov test.
I do not understand how these tests could help me with a goodnessoffit test for the Von MisesFisher distribution, but probably is my limit.
Could anyone being of any help?
Regards,
Sergio
Marc (view profile)
sergio  did you see the pdf with descriptions? (http://www.jstatsoft.org/v31/i10/paper)
You probably want either the ktest of the kuipertest.
sergio (view profile)
Hi guys, I'm new to circular statistics and I've downloaded this package.
Given some vectors, I'd like to test if they are distributed following a Von MisesFisher distribution.
Do you know what instructions of the package I should use?
Can you help?
Aviram Gelblum (view profile)
After some testing I figured the previous bug has to do with recurrence of unique values in the data. I took care of it by using
alpha=alpha+0.00001*(1:numel(alpha)), but this is obviously a workaround which isn't satisfactory for a selfrespecting algorithm.
At any rate, I forgot to mention how great this toolbox is. It has been of great help, and saved me a lot of time and work.
Aviram Gelblum (view profile)
Hi, I'm getting wrong clustering using circ_clust.
for example, if I give as an input circ_clust([1 1 1 1 3.5 4 5.5 0.5],2)
I get
ans =
1
2
1
1
1
1
1
1
sometimes the clustering does work, but I don't know why it does/doesn't...
I'm using Matlab 2012b...
Dylan Muir (view profile)
Great toolbox. I was wondering if it is possible to have more accurate pvalue estimates in the Kuiper test?
Francesco Montorsi (view profile)
Hi,
great toolbox, thanks.
By the way, I agree with Allan's comment (see below) that regarding the Von Mises distribution, it may be useful to have an implementation with higher numerical stability. In particular, I added this trivial function, which returns the logpdf of the Von Mises distrib:
function [p alpha] = circ_vm_logpdf(alpha, thetahat, kappa)
% if no angles are supplied, 100 evenly spaced points around the circle are
% chosen
if nargin < 1  isempty(alpha)
alpha = linspace(0, 2*pi, 101)';
alpha = alpha(1:end1);
end
if nargin < 3
kappa = 1;
end
if nargin < 2
thetahat = 0;
end
alpha = alpha(:);
% evaluate pdf
C = log( 2*pi*besseli(0,kappa) );
p = C + kappa*cos(alphathetahat);
Thanks to the greater numerical stability logpdfs are often used in place of pdfs, so this little function may be of help to others...
Allan (view profile)
Hi there, great toolbox. I propose a change to avoid numerical instability in circ_vmpdf.m.
Current code to evaluate the pdf:
C = 1/(2*pi*besseli(0,kappa));
p = C * exp(kappa*cos(alphathetahat));
Proposed replacement code:
C = log(1)log(2*pi*besseli(0,kappa,1))+(kappa*cos(alphathetahat))kappa;
p = exp(C);
Examples:
circ_vmpdf(0,0,1000)
Old code result: NaN
New code result: 12.6141
Philipp Berens (view profile)
Ryan, the average is in the dot product w'*exp(...) which in the simplest case is a vector of ones  so this is the sum operation. exp(i*angle) decomposes the angle into its sine and cosine components. Finally, angle is atan2. Compare the results of your and my code  they should be identical with my code likely running a bit fast due to matrix style computations.
Bst
Philipp
Ryan (view profile)
I haven't run through this toolbox yet, so I apologize if I am missing something with this question (I just glanced through the source code because I am interested in directional stats).
When you calculate the mean, the formula you use is:
% compute weighted sum of cos and sin of angles
r = w'*exp(1i*alpha);
% obtain mean by
mu = angle(r);
Now, correct me if I'm wrong, but this doesn't seem to calculate the average at all? It seems to me that here we are inputting a data array into the angle command, which will output the phase angle of each element of that array, not a singular mean.
Wouldn't a better way of calculating the average be to use atan2? Something like:
for i = 1:w
S(i) = sin(alpha(i));
C(i) = cos(alpha(i));
end
X = sum(S)*(1/w);
Y = sum(C)*(1/w);
mu = atan2(X,Y);
Marnix Maas (view profile)
Thanks for the great toolbox! I have a question: I have a set of directional stochastic variables that are mutually correlated. I have used circ_corrcc to construct a correlation matrix for these variables, but I’m also interested in their covariance matrix. There does not appear to be a function for this in the current toolbox.
Not having any previous experience with circular statistics, I’m wondering if it makes sense to construct a covariance matrix by denormalizing the correlation matrix, multiplying each element by the two corresponding circular standard deviations? Perhaps a covariance matrix could be a useful addition to the toolbox.
Thanks,
Marnix
Philipp Berens (view profile)
Hi Francesco, if you have orientations, multiply all orientations by 2 to obtain directions. If you want to obtain the mean resultant vector, devide its orientation by 2 again.
Francesco (view profile)
excuse my previous post! I just realize what paxial truly meant.
For further reference this will solve the previously cited problem
%% uniform distribution test
% in the interval [0, 180)
y180 = circ_axial(circ_ang2rad(0 + 179*rand(4000,1)),2);
p180 = circ_otest(y180)
% in the interval [0 360)
y360 = deg2rad(0 + 359*rand(4000,1));
p360 = circ_otest(y360)
Francesco (view profile)
excuse my previous post! I just realize what paxial truly meant.
For further reference this will solve the previously cited problem
%% uniform distribution test
% in the interval [0, 180)
y180 = circ_axial(circ_ang2rad(0 + 179*rand(4000,1)),2);
p180 = circ_otest(y180)
% in the interval [0 360)
y360 = deg2rad(0 + 359*rand(4000,1));
p360 = circ_otest(y360)
Francesco (view profile)
I am testing the toolbox out with not much of a prior knowledge on the subject. It seems a really good piece of software and it's helping me out grasping some of the theory.
I have a question: if I am dealing with orientations [0, 180) degrees more than directions [0 360), is there a proper way to transform may data prior to using the function in the toolbox?
For example, if I am trying to test for circular uniformity with a population that is uniformly distributed in [0 180)  which I'd like to have a p>0.05  I obtained a very small value, which is consistent with the test looking over the full interval.
Suggestions? Thanks
Francesco
Example code 
y180 = circ_ang2rad(0 + 179*rand(4000,1));
p180 = circ_otest(y180)
% in the interval [0 360)
y360 = deg2rad(0 + 359*rand(4000,1));
p360 = circ_otest(y360)
Jer Walley (view profile)
Great toolbox! Exactly what I needed. However, my data has many NaN's  do you have a way to work around data with gaps?
Matt Davis (view profile)
This is a great toolbox  very helpful. A few bug reports:
1. formatSubPlot calls "parseVarArgs" that's not standard matlab, or part of this toolbox. Could you add a pointer to where to download this.
2. In Example 2 the descriptive stats cell needs updating to respect the matrix style computations. So, line 67 should read:
stats(i,1) = circ_mean(ori,spk,2);
and similar for all the other lines of code.
Thanks for supporting this toolbox.
Natalia (view profile)
I confirm Dillon's report on circ_wwtest bug.
Dillon (view profile)
Great toolbox but I think there is an error in the logic used at circ_wwtest.m > checkAssumption() lines 107115.
The code currently reads:
if n > 10 && rw<.45
warning('Test not applicable. Average resultant vector length < 0.45.') %#ok<WNTAG>
elseif n > 6 && rw<.5
warning('Test not applicable. Average number of samples per population < 11 and average resultant vector length < 0.5.') %#ok<WNTAG>
elseif n >=5 && rw<.55
warning('Test not applicable. Average number of samples per population < 7 and average resultant vector length < 0.55.') %#ok<WNTAG>
elseif n < 5
warning('Test not applicable. Average number of samples per population < 5.') %#ok<WNTAG>
end
Notice that the if/else statements do not match the warning text. Particularly when n>5 the user will always be warned when the resultant vector, rw<0.55 which is not captured by the warning. The corrected if/else statements are as follows:
if n >= 11 && rw<.45
warning('Test not applicable. Average resultant vector length < 0.45.') %#ok<WNTAG>
elseif n<11 && n >= 7 && rw<.5
warning('Test not applicable. Average number of samples per population < 11 and average resultant vector length < 0.5.') %#ok<WNTAG>
elseif n<7 && n >=5 && rw<.55
warning('Test not applicable. Average number of samples per population < 7 and average resultant vector length < 0.55.') %#ok<WNTAG>
elseif n < 5
warning('Test not applicable. Average number of samples per population < 5.') %#ok<WNTAG>
end
I've assumed that the warning statements are correct but if the if/else statements are correct it would be more compact to warn the user under only 2 conditions: n<5 and rw<0.55.
Thanks again for the very useful toolbox.
Philipp Berens (view profile)
Thanks for the recent feedback and bugreports. I was away for a while and will start taking care of them soon.
Omar Mian (view profile)
Suggestion for addition:
Parametric and nonparametric paired sample tests, Zar (2010) Biostatistical Analysis, sections 27.13 and 27.14
Pablo Martinez (view profile)
Sorry wrong line number. The error is at line 169!
in the function circ_hktest.m
pI = 1  chi2pdf(chiI, df_i);
It should be
pI = 1  chi2cdf(chiI, df_i);
Great toolbox.
Pablo Martinez (view profile)
Great toolbox!
Pablo Martinez (view profile)
I found an error in the function circ_hktest.m at line 160
pI = 1  chi2pdf(chiI, df_i);
It should be
pI = 1  chi2cdf(chiI, df_i);
Natalia (view profile)
Thank you very much for such a useful toolbox. Now, I have a question related to circ_ktest (twosmple test to compare concentration). The F statistic is defined only in case of rbar>.7, Mardia (pag 133, 1999) compute F in the case where resultant vector length is <0.45 :
n1 = length(alpha1);
n2 = length(alpha2);
R1avg=circ_r(alpha1);
R2avg = circ_r(alpha2);
R1 = n1*circ_r(alpha1);
R2 = n2*circ_r(alpha2);
%make sure that rbar > .7
rbar = (R1+R2)/(n1+n2);
if rbar > .7
f = ((n21)*(n1R1))/((n11)*(n2R2));
elseif rbar< .45 %taken from Mardia 1999 p.133 (Baschelet report: Mardia 1972 pag 161)
g11= asin(2*sqrt(3/8)*(R1avg));
g12= asin(2*sqrt(3/8)*(R2avg));
f= (2/sqrt(3))*((g11g12)/(1/(n14)+ 1/(n24)).^(1/2));
But here Sample 1 and Sample 2 define the sign of F... and so S1 and S2 will be defined depending on Ravg value being S1>S2 for computation of F. Is this right?
Thank you!
natalia
Luke (view profile)
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]
Balazs Barkoczi (view profile)
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!
FuhCherng (view profile)
@Christopher: Thank you so much for your kindness and help. I really appreciate it.
Christopher (view profile)
@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.
FuhCherng (view profile)
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
Philipp Berens (view profile)
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 multisample tests assume independent samples. I don't know about repeatedmeasures 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.
Omar Mian (view profile)
Can any of the multisample tests in this toolbox be used with repeated measures data or do they all assume independent samples?
Omar Mian (view profile)
Very useful toolbox. Option to ignore NaNs in the calculations would make it even better.
Vlad Atanasiu (view profile)
Good toolbox! I added a function for kernel smoothing density estimate for circular data here: http://www.mathworks.com/matlabcentral/fileexchange/32614kernelsmoothingdensityestimateforcirculardata .
Heida Maria Sigurdardottir (view profile)
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 nonparametric 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 pvalue correction?
Christopher (view profile)
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.
Christopher (view profile)
Sorry, the error is in circ_kuipertest.m, not circ_kuiper.m as previously written!
Quang Nguyen Dinh (view profile)
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 misunderstand 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.
Bart Geurten (view profile)
Needed a circular statistic means and got what I needed... thanks!
Philipp Berens (view profile)
Hi Christopher, thanks for your feedback. I will update circ_plot with the next upload.
Omar Mian (view profile)
Christopher (view profile)
Another issue with circ_plot.m:
A typo in line 121, should read
s = varargin{3};
(instead of vargin{1})
Christopher (view profile)
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 (view profile)
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).
Carlos (view profile)
Nice work
Kourosh (view profile)
Excellent toolbox!!! Thanks Philipp!
Philipp Berens (view profile)
I fixed the bug in circ_clust.
Aisling (view profile)
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 (view profile)
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:end1).', mu, kappa);
p = p * diff(points(1:2));
to obtain approximate probabilities.
Cesare (view profile)
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:end1).', 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
Cesare (view profile)
also with the 2010b.
Thanks,
Cesare
Philipp Berens (view profile)
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 'Matlabstyle'. 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.
Tomo (view profile)
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.
Cesare (view profile)
something's wrong in the new circ_vmpdf....can't replicate results from previous release
chairmanK (view profile)
Functions do not gracefully handle Ndimensional arrays for powerful MATLABstyle 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!
Philipp Berens (view profile)
Issue is fixed in the upload of 11/5/09.
Omar Mian (view profile)
Encountering same issue as Shiquan
Shiquan Wang (view profile)
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'.
Shiquan Wang (view profile)
This is great work! Thanks.
Philipp Berens (view profile)
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
Flurin Honegger (view profile)
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/(11/(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 (view profile)
Both definitions are around... I will optionally add computing both with the next update.
Richard Heitz (view profile)
Sorry, but I think circ_var returns s = (1r) when it should be s = 2*(1r).
for this reason, circ_std and circ_var will not agree
Richard Courtemanche (view profile)
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!
Tal Krasovsky (view profile)
Excellent toolbox, helped me a lot. Greatly appreciated!
Richard Heitz (view profile)
this is a great toolbox
Lorenzo Guerrasio (view profile)
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
Lorenzo Guerrasio (view profile)
I wish I saw it before :)
Very nicely done.
I take that back  guilty of confusing 'r' and 'R'
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;
Would include skewness and kurtosis!
Exactly what I needed
thanks.