spm_Bcdf.m 3.0 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091
  1. function F = spm_Bcdf(x,v,w)
  2. % Inverse Cumulative Distribution Function (CDF) of Beta distribution
  3. % FORMAT F = spm_Bcdf(x,v,w)
  4. %
  5. % x - Beta variates (Beta has range [0,1])
  6. % v - Shape parameter (v>0)
  7. % w - Shape parameter (w>0)
  8. % F - CDF of Beta distribution with shape parameters [v,w] at points x
  9. %__________________________________________________________________________
  10. %
  11. % spm_Bcdf implements the Cumulative Distribution Function for Beta
  12. % distributions.
  13. %
  14. % Definition:
  15. %--------------------------------------------------------------------------
  16. % The Beta distribution has two shape parameters, v and w, and is
  17. % defined for v>0 & w>0 and for x in [0,1] (See Evans et al., Ch5).
  18. % The Cumulative Distribution Function (CDF) F(x) is the probability
  19. % that a realisation of a Beta random variable X has value less than
  20. % x. F(x)=Pr{X<x}: This function is usually known as the incomplete Beta
  21. % function. See Abramowitz & Stegun, 26.5; Press et al., Sec6.4 for
  22. % definitions of the incomplete beta function.
  23. %
  24. % Variate relationships:
  25. %--------------------------------------------------------------------------
  26. % Many: See Evans et al., Ch5
  27. %
  28. % Algorithm:
  29. %--------------------------------------------------------------------------
  30. % Using MATLAB's implementation of the incomplete beta finction (betainc).
  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. % Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging
  47. % Andrew Holmes
  48. % $Id: spm_Bcdf.m 4182 2011-02-01 12:29:09Z guillaume $
  49. %-Format arguments, note & check sizes
  50. %--------------------------------------------------------------------------
  51. if nargin<3, error('Insufficient arguments'), end
  52. ad = [ndims(x);ndims(v);ndims(w)];
  53. rd = max(ad);
  54. as = [[size(x),ones(1,rd-ad(1))];...
  55. [size(v),ones(1,rd-ad(2))];...
  56. [size(w),ones(1,rd-ad(3))]];
  57. rs = max(as);
  58. xa = prod(as,2)>1;
  59. if sum(xa)>1 && any(any(diff(as(xa,:)),1))
  60. error('non-scalar args must match in size');
  61. end
  62. %-Computation
  63. %--------------------------------------------------------------------------
  64. %-Initialise result to zeros
  65. F = zeros(rs);
  66. %-Only defined for x in [0,1] & strictly positive v & w.
  67. % Return NaN if undefined.
  68. md = ( x>=0 & x<=1 & v>0 & w>0 );
  69. if any(~md(:))
  70. F(~md) = NaN;
  71. warning('Returning NaN for out of range arguments');
  72. end
  73. %-Special cases: F=1 when x=1
  74. F(md & x==1) = 1;
  75. %-Non-zero where defined & x>0, avoid special cases
  76. Q = find( md & x>0 & x<1 );
  77. if isempty(Q), return, end
  78. if xa(1), Qx=Q; else Qx=1; end
  79. if xa(2), Qv=Q; else Qv=1; end
  80. if xa(3), Qw=Q; else Qw=1; end
  81. %-Compute
  82. F(Q) = betainc(x(Qx),v(Qv),w(Qw));