spm_Ncdf.m 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293
  1. function F = spm_Ncdf(x,u,v)
  2. % Cumulative Distribution Function (CDF) for univariate Normal distributions
  3. % FORMAT F = spm_Ncdf(x,u,v)
  4. %
  5. % x - ordinates
  6. % u - mean [Defaults to 0]
  7. % v - variance (v>0) [Defaults to 1]
  8. % F - pdf of N(u,v) at x (Lower tail probability)
  9. %__________________________________________________________________________
  10. %
  11. % spm_Ncdf implements the Cumulative Distribution Function (CDF) for
  12. % the Normal (Gaussian) family of distributions.
  13. %
  14. % Definition:
  15. %--------------------------------------------------------------------------
  16. % The CDF F(x) of a Normal distribution with mean u and variance v is
  17. % the probability that a random realisation X from this distribution
  18. % will be less than x. F(x)=Pr(X<=x) for X~N(u,v). See Evans et al.,
  19. % Ch29 for further definitions and variate relationships.
  20. %
  21. % If X~N(u,v), then Z=(Z-u)/sqrt(v) has a standard normal distribution,
  22. % Z~N(0,1). The CDF of the standard normal distribution is known as \Phi(z).
  23. %
  24. % (KWorsley) For extreme variates with abs(z)>6 where z=(x-u)/sqrt(v), the
  25. % approximation \Phi(z) \approx exp(-z^2/2)/(z*sqrt(2*pi)) may be useful.
  26. %
  27. % Algorithm:
  28. %--------------------------------------------------------------------------
  29. % The CDF for a standard N(0,1) Normal distribution, \Phi(z), is
  30. % related to the error function by: (Abramowitz & Stegun, 26.2.29)
  31. %
  32. % \Phi(z) = 0.5 + erf(z/sqrt(2))/2
  33. %
  34. % MATLAB's implementation of the error function is used for computation.
  35. %
  36. % References:
  37. %--------------------------------------------------------------------------
  38. % Evans M, Hastings N, Peacock B (1993)
  39. % "Statistical Distributions"
  40. % 2nd Ed. Wiley, New York
  41. %
  42. % Abramowitz M, Stegun IA, (1964)
  43. % "Handbook of Mathematical Functions"
  44. % US Government Printing Office
  45. %
  46. % Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
  47. % "Numerical Recipes in C"
  48. % Cambridge
  49. %
  50. %__________________________________________________________________________
  51. % Copyright (C) 1995-2011 Wellcome Trust Centre for Neuroimaging
  52. % Andrew Holmes
  53. % $Id: spm_Ncdf.m 4182 2011-02-01 12:29:09Z guillaume $
  54. %-Format arguments, note & check sizes
  55. %--------------------------------------------------------------------------
  56. if nargin<3, v=1; end
  57. if nargin<2, u=0; end
  58. if nargin<1, F=[]; return, end
  59. ad = [ndims(x);ndims(u);ndims(v)];
  60. rd = max(ad);
  61. as = [[size(x),ones(1,rd-ad(1))];...
  62. [size(u),ones(1,rd-ad(2))];...
  63. [size(v),ones(1,rd-ad(3))]];
  64. rs = max(as);
  65. xa = prod(as,2)>1;
  66. if sum(xa)>1 && any(any(diff(as(xa,:)),1))
  67. error('non-scalar args must match in size');
  68. end
  69. %-Computation
  70. %--------------------------------------------------------------------------
  71. %-Initialise result to zeros
  72. F = zeros(rs);
  73. %-Only defined for strictly positive variance v. Return NaN if undefined.
  74. md = ( ones(size(x)) & ones(size(u)) & v>0 );
  75. if any(~md(:))
  76. F(~md) = NaN;
  77. warning('SPM:negativeVariance','Returning NaN for out of range arguments.');
  78. end
  79. %-Non-zero where defined
  80. Q = find( md );
  81. if isempty(Q), return, end
  82. if xa(1), Qx=Q; else Qx=1; end
  83. if xa(2), Qu=Q; else Qu=1; end
  84. if xa(3), Qv=Q; else Qv=1; end
  85. %-Compute
  86. F(Q) = 0.5 + 0.5*erf((x(Qx)-u(Qu))./sqrt(2*v(Qv)));