bootmedian.m 1.8 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586
  1. function varargout=bootmedian(varargin)
  2. % [p,ci]=bootmedian(A,B,['boots'])
  3. % When passed in a single vector, tests whether the median is different from 0
  4. % When passed in two vectors does a permutation test to see whether the medians
  5. % are significantly different
  6. % Optional inputs:
  7. % 'boots' the number of shuffles to perform (2000)
  8. if nargin==0
  9. help('stats.bootmedian')
  10. return
  11. end
  12. if nargin>1 && isnumeric(varargin{2})
  13. A=varargin{1};
  14. B=varargin{2};
  15. one_dist=false;
  16. if nargin>2
  17. varargin=varargin(3:end);
  18. else
  19. varargin={};
  20. end
  21. else
  22. A=varargin{1};
  23. one_dist=true;
  24. if nargin>1
  25. varargin=varargin(2:end);
  26. else
  27. varargin={};
  28. end
  29. end
  30. boots=2000;
  31. utils.overridedefaults(who,varargin);
  32. if one_dist
  33. % assume test whether mean of the population is differenct from zero.
  34. dist=A;
  35. if isvector(dist)
  36. dist=dist(:);
  37. end
  38. n=size(dist,1);
  39. [B]=bootstrp(boots, @nanmedian, dist);
  40. ps=[0:0.01:100];
  41. sd_ps=prctile(B,ps);
  42. sd_p=stats.get_p(0,B);
  43. varargout{1}=sd_p;
  44. varargout{2}=prctile(B,50);
  45. varargout{3}=B;
  46. elseif nargin==2
  47. sA=size(A,1);
  48. sB=size(B,1);
  49. sd=nanmedian(A)-nanmedian(B);
  50. if min(sA,sB)<=6
  51. warning('Not meaningful to compute 2000 bootstraps when n<7');
  52. boots=factorial(min(sA,sB));
  53. end
  54. ALL_DATA=[A;B];
  55. boot_score=zeros(boots,size(ALL_DATA,2));
  56. for bx=1:boots
  57. shuff_d=ALL_DATA(randperm(sA+sB),:);
  58. A=shuff_d(1:sA,:);
  59. B=shuff_d(sA+1:end,:);
  60. boot_score(bx,:)=nanmedian(A)-nanmedian(B);
  61. end
  62. sd_p=stats.get_p(sd, boot_score);
  63. end
  64. varargout{1}=sd_p;
  65. varargout{2}=prctile(B,[2.5 97.5]);
  66. varargout{3}=prctile(B,[50]);