circ_rtest.m 1.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475
  1. function [pval z] = circ_rtest(alpha, w, d)
  2. %
  3. % [pval, z] = circ_rtest(alpha,w)
  4. % Computes Rayleigh test for non-uniformity of circular data.
  5. % H0: the population is uniformly distributed around the circle
  6. % HA: the populatoin is not distributed uniformly around the circle
  7. % Assumption: the distribution has maximally one mode and the data is
  8. % sampled from a von Mises distribution!
  9. %
  10. % Input:
  11. % alpha sample of angles in radians
  12. % [w number of incidences in case of binned angle data]
  13. % [d spacing of bin centers for binned data, if supplied
  14. % correction factor is used to correct for bias in
  15. % estimation of r, in radians (!)]
  16. %
  17. % Output:
  18. % pval p-value of Rayleigh's test
  19. % z value of the z-statistic
  20. %
  21. % PHB 7/6/2008
  22. %
  23. % References:
  24. % Statistical analysis of circular data, N. I. Fisher
  25. % Topics in circular statistics, S. R. Jammalamadaka et al.
  26. % Biostatistical Analysis, J. H. Zar
  27. %
  28. % Circular Statistics Toolbox for Matlab
  29. % By Philipp Berens, 2009
  30. % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
  31. if size(alpha,2) > size(alpha,1)
  32. alpha = alpha';
  33. end
  34. if nargin < 2
  35. r = circ_r(alpha);
  36. n = length(alpha);
  37. else
  38. if length(alpha)~=length(w)
  39. error('Input dimensions do not match.')
  40. end
  41. if nargin < 3
  42. d = 0;
  43. end
  44. r = circ_r(alpha,w(:),d);
  45. n = sum(w);
  46. end
  47. % compute Rayleigh's R (equ. 27.1)
  48. R = n*r;
  49. % compute Rayleigh's z (equ. 27.2)
  50. z = R^2 / n;
  51. % compute p value using approxation in Zar, p. 617
  52. pval = exp(sqrt(1+4*n+4*(n^2-R^2))-(1+2*n));
  53. % outdated version:
  54. % compute the p value using an approximation from Fisher, p. 70
  55. % pval = exp(-z);
  56. % if n < 50
  57. % pval = pval * (1 + (2*z - z^2) / (4*n) - ...
  58. % (24*z - 132*z^2 + 76*z^3 - 9*z^4) / (288*n^2));
  59. % end