circ_confmean.m 2.0 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
  1. function t = circ_confmean(alpha, xi, w, d, dim)
  2. %
  3. % t = circ_mean(alpha, xi, w, d, dim)
  4. % Computes the confidence limits on the mean for circular data.
  5. %
  6. % Input:
  7. % alpha sample of angles in radians
  8. % [xi (1-xi)-confidence limits are computed, default 0.05]
  9. % [w number of incidences in case of binned angle data]
  10. % [d spacing of bin centers for binned data, if supplied
  11. % correction factor is used to correct for bias in
  12. % estimation of r, in radians (!)]
  13. % [dim compute along this dimension, default is 1]
  14. %
  15. % Output:
  16. % t mean +- d yields upper/lower (1-xi)% confidence limit
  17. %
  18. % PHB 7/6/2008
  19. %
  20. % References:
  21. % Statistical analysis of circular data, N. I. Fisher
  22. % Topics in circular statistics, S. R. Jammalamadaka et al.
  23. % Biostatistical Analysis, J. H. Zar
  24. %
  25. % Circular Statistics Toolbox for Matlab
  26. % By Philipp Berens, 2009
  27. % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
  28. if nargin < 5
  29. dim = 1;
  30. end
  31. if nargin < 4 || isempty(d)
  32. % per default do not apply correct for binned data
  33. d = 0;
  34. end
  35. if nargin < 3 || isempty(w)
  36. % if no specific weighting has been specified
  37. % assume no binning has taken place
  38. w = ones(size(alpha));
  39. else
  40. if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
  41. error('Input dimensions do not match');
  42. end
  43. end
  44. % set confidence limit size to default
  45. if nargin < 2 || isempty(xi)
  46. xi = 0.05;
  47. end
  48. % compute ingredients for conf. lim.
  49. r = circ_r(alpha,w,d,dim);
  50. n = sum(w,dim);
  51. R = n.*r;
  52. c2 = chi2inv((1-xi),1);
  53. % check for resultant vector length and select appropriate formula
  54. t = zeros(size(r));
  55. for i = 1:numel(r)
  56. if r(i) < .9 && r(i) > sqrt(c2/2/n(i))
  57. t(i) = sqrt((2*n(i)*(2*R(i)^2-n(i)*c2))/(4*n(i)-c2)); % equ. 26.24
  58. elseif r(i) >= .9
  59. t(i) = sqrt(n(i)^2-(n(i)^2-R(i)^2)*exp(c2/n(i))); % equ. 26.25
  60. else
  61. t(i) = NaN;
  62. warning('Requirements for confidence levels not met.');
  63. end
  64. end
  65. % apply final transform
  66. t = acos(t./R);