spm_invIcdf.m 3.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495
  1. function r = spm_invIcdf(F,n,p)
  2. % Inverse Cumulative Distribution Function (CDF) of Binomial distribution
  3. % FORMAT r = spm_invIcdf(F,n,p)
  4. %
  5. % F - CDF (lower tail p-value)
  6. % n - Binomial n
  7. % p - Binomial p [Defaults to 0.5]
  8. % r - ordinate
  9. %__________________________________________________________________________
  10. %
  11. % spm_invIcdf returns the inverse Cumulative Distribution Function for
  12. % the Binomial family of distributions.
  13. %
  14. % Definition:
  15. %--------------------------------------------------------------------------
  16. % The Bin(n,p) distribution is the distribution of the number of
  17. % successes from n identical independent Bernoulli trials each with
  18. % success probability p. If random variable R is the number of
  19. % successes from such a set of Bernoulli trials, then the CDF F(r) is
  20. % the probability of r or less sucesses.
  21. %
  22. % The Binomial distribution is discrete, defined for p in [0,1] and r
  23. % in {0,1,...,n}, so F(r) is a discrete function. This inverse CDF
  24. % function returns the smallest Whole r such that the F(r) equals or
  25. % exceeds the given CDF probability F. I.e. F(r) is treated as a step
  26. % function.
  27. %
  28. % Algorithm:
  29. %--------------------------------------------------------------------------
  30. % r is found by direct summation of the Binomial PDFs until F is exceeded.
  31. %
  32. % References:
  33. %--------------------------------------------------------------------------
  34. % Evans M, Hastings N, Peacock B (1993)
  35. % "Statistical Distributions"
  36. % 2nd Ed. Wiley, New York
  37. %
  38. % Abramowitz M, Stegun IA, (1964)
  39. % "Handbook of Mathematical Functions"
  40. % US Government Printing Office
  41. %
  42. % Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
  43. % "Numerical Recipes in C"
  44. % Cambridge
  45. %
  46. %__________________________________________________________________________
  47. % Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging
  48. % Andrew Holmes
  49. % $Id: spm_invIcdf.m 4182 2011-02-01 12:29:09Z guillaume $
  50. %-Format arguments, note & check sizes
  51. %--------------------------------------------------------------------------
  52. if nargin<3, p=0.5; end
  53. if nargin<2, error('Insufficient arguments'), end
  54. ad = [ndims(F);ndims(n);ndims(p)];
  55. rd = max(ad);
  56. as = [[size(F),ones(1,rd-ad(1))];...
  57. [size(n),ones(1,rd-ad(2))];...
  58. [size(p),ones(1,rd-ad(3))]];
  59. rs = max(as);
  60. xa = prod(as,2)>1;
  61. if sum(xa)>1 && any(any(diff(as(xa,:)),1))
  62. error('non-scalar args must match in size');
  63. end
  64. %-Computation
  65. %--------------------------------------------------------------------------
  66. %-Initialise result to zeros
  67. r = zeros(rs);
  68. %-Only defined for whole n, p&F in [0,1]. Return NaN if undefined.
  69. md = ( F>=0 & F<=1 & n==floor(n) & n>=0 & p>=0 & p<=1 );
  70. if any(~md(:))
  71. r(~md) = NaN;
  72. warning('Returning NaN for out of range arguments');
  73. end
  74. %-Non-zero only where defined, & F>0
  75. Q = find( md & F>0 );
  76. if isempty(Q), return, end
  77. if xa(1), QF=Q; else QF=1; end
  78. if xa(2), Qn=Q; else Qn=1; end
  79. if xa(3), Qp=Q; else Qp=1; end
  80. %-Compute by directly summing Bin PDF's for successive r & comparing with F
  81. tr = 0;
  82. Ftr = spm_Ipdf(tr,n(Qn),p(Qp));
  83. while any(F(QF)>Ftr) && any(n(Qn)>tr)
  84. tr = tr+1;
  85. i = find(Ftr<F(QF));
  86. r(Q(i)) = r(Q(i)) + 1;
  87. Ftr = Ftr + spm_Ipdf(tr,n(Qn),p(Qp));
  88. end