octave_packages/statistics-1.1.3/anderson_darling_cdf.m

   1 ## Author: Paul Kienzle <pkienzle@users.sf.net>
   2 ## This program is granted to the public domain.
   3
   4 ## -*- texinfo -*-
   5 ## @deftypefn {Function File} @var{p} = anderson_darling_cdf (@var{A}, @var{n})
   6 ##
   7 ## Return the CDF for the given Anderson-Darling coefficient @var{A}
   8 ## computed from @var{n} values sampled from a distribution. For a
   9 ## vector of random variables @var{x} of length @var{n}, compute the CDF
  10 ## of the values from the distribution from which they are drawn.
  11 ## You can uses these values to compute @var{A} as follows:
  12 ##
  13 ## @example
  14 ## @var{A} = -@var{n} - sum( (2*i-1) .* (log(@var{x}) + log(1 - @var{x}(@var{n}:-1:1,:))) )/@var{n};
  15 ## @end example
  16 ##
  17 ## From the value @var{A}, @code{anderson_darling_cdf} returns the probability
  18 ## that @var{A} could be returned from a set of samples.
  19 ##
  20 ## The algorithm given in [1] claims to be an approximation for the
  21 ## Anderson-Darling CDF accurate to 6 decimal points.
  22 ##
  23 ## Demonstrate using:
  24 ##
  25 ## @example
  26 ## n = 300; reps = 10000;
  27 ## z = randn(n, reps);
  28 ## x = sort ((1 + erf (z/sqrt (2)))/2);
  29 ## i = [1:n]' * ones (1, size (x, 2));
  30 ## A = -n - sum ((2*i-1) .* (log (x) + log (1 - x (n:-1:1, :))))/n;
  31 ## p = anderson_darling_cdf (A, n);
  32 ## hist (100 * p, [1:100] - 0.5);
  33 ## @end example
  34 ##
  35 ## You will see that the histogram is basically flat, which is to
  36 ## say that the probabilities returned by the Anderson-Darling CDF
  37 ## are distributed uniformly.
  38 ##
  39 ## You can easily determine the extreme values of @var{p}:
  40 ##
  41 ## @example
  42 ## [junk, idx] = sort (p);
  43 ## @end example
  44 ##
  45 ## The histograms of various @var{p} aren't  very informative:
  46 ##
  47 ## @example
  48 ## histfit (z (:, idx (1)), linspace (-3, 3, 15));
  49 ## histfit (z (:, idx (end/2)), linspace (-3, 3, 15));
  50 ## histfit (z (:, idx (end)), linspace (-3, 3, 15));
  51 ## @end example
  52 ##
  53 ## More telling is the qqplot:
  54 ##
  55 ## @example
  56 ## qqplot (z (:, idx (1))); hold on; plot ([-3, 3], [-3, 3], ';;'); hold off;
  57 ## qqplot (z (:, idx (end/2))); hold on; plot ([-3, 3], [-3, 3], ';;'); hold off;
  58 ## qqplot (z (:, idx (end))); hold on; plot ([-3, 3], [-3, 3], ';;'); hold off;
  59 ## @end example
  60 ##
  61 ## Try a similarly analysis for @var{z} uniform:
  62 ##
  63 ## @example
  64 ## z = rand (n, reps); x = sort(z);
  65 ## @end example
  66 ##
  67 ## and for @var{z} exponential:
  68 ##
  69 ## @example
  70 ## z = rande (n, reps); x = sort (1 - exp (-z));
  71 ## @end example
  72 ##
  73 ## [1] Marsaglia, G; Marsaglia JCW; (2004) "Evaluating the Anderson Darling
  74 ## distribution", Journal of Statistical Software, 9(2).
  75 ##
  76 ## @seealso{anderson_darling_test}
  77 ## @end deftypefn
  78
  79 function y = anderson_darling_cdf(z,n)
  80   y = ADinf(z);
  81   y += ADerrfix(y,n);
  82 end
  83
  84 function y = ADinf(z)
  85   y = zeros(size(z));
  86
  87   idx = (z < 2);
  88   if any(idx(:))
  89     p = [.00168691, -.0116720, .0347962, -.0649821, .247105, 2.00012];
  90     z1 = z(idx);
  91     y(idx) = exp(-1.2337141./z1)./sqrt(z1).*polyval(p,z1);
  92   end
  93
  94   idx = (z >= 2);
  95   if any(idx(:))
  96     p = [-.0003146, +.008056, -.082433, +.43424, -2.30695, 1.0776];
  97     y(idx) = exp(-exp(polyval(p,z(idx))));
  98   end
  99 end
 100
 101 function y = ADerrfix(x,n)
 102   if isscalar(n), n = n*ones(size(x));
 103   elseif isscalar(x), x = x*ones(size(n));
 104   end
 105   y = zeros(size(x));
 106   c = .01265 + .1757./n;
 107
 108   idx = (x >= 0.8);
 109   if any(idx(:))
 110     p = [255.7844, -1116.360, 1950.646, -1705.091, 745.2337, -130.2137];
 111     g3 = polyval(p,x(idx));
 112     y(idx) = g3./n(idx);
 113   end
 114
 115   idx = (x < 0.8 & x > c);
 116   if any(idx(:))
 117     p = [1.91864, -8.259, 14.458, -14.6538, 6.54034, -.00022633];
 118     n1 = 1./n(idx);
 119     c1 = c(idx);
 120     g2 = polyval(p,(x(idx)-c1)./(.8-c1));
 121     y(idx) = (.04213 + .01365*n1).*n1 .* g2;
 122   end
 123
 124   idx = (x <= c);
 125   if any(idx(:))
 126     x1 = x(idx)./c(idx);
 127     n1 = 1./n(idx);
 128     g1 = sqrt(x1).*(1-x1).*(49*x1-102);
 129     y(idx) = ((.0037*n1+.00078).*n1+.00006).*n1 .* g1;
 130   end
 131 end