1 ## Author: Paul Kienzle <pkienzle@users.sf.net>
2 ## This program is granted to the public domain.
5 ## @deftypefn {Function File} @var{p} = anderson_darling_cdf (@var{A}, @var{n})
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:
14 ## @var{A} = -@var{n} - sum( (2*i-1) .* (log(@var{x}) + log(1 - @var{x}(@var{n}:-1:1,:))) )/@var{n};
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.
20 ## The algorithm given in [1] claims to be an approximation for the
21 ## Anderson-Darling CDF accurate to 6 decimal points.
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);
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.
39 ## You can easily determine the extreme values of @var{p}:
42 ## [junk, idx] = sort (p);
45 ## The histograms of various @var{p} aren't very informative:
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));
53 ## More telling is the qqplot:
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;
61 ## Try a similarly analysis for @var{z} uniform:
64 ## z = rand (n, reps); x = sort(z);
67 ## and for @var{z} exponential:
70 ## z = rande (n, reps); x = sort (1 - exp (-z));
73 ## [1] Marsaglia, G; Marsaglia JCW; (2004) "Evaluating the Anderson Darling
74 ## distribution", Journal of Statistical Software, 9(2).
76 ## @seealso{anderson_darling_test}
79 function y = anderson_darling_cdf(z,n)
89 p = [.00168691, -.0116720, .0347962, -.0649821, .247105, 2.00012];
91 y(idx) = exp(-1.2337141./z1)./sqrt(z1).*polyval(p,z1);
96 p = [-.0003146, +.008056, -.082433, +.43424, -2.30695, 1.0776];
97 y(idx) = exp(-exp(polyval(p,z(idx))));
101 function y = ADerrfix(x,n)
102 if isscalar(n), n = n*ones(size(x));
103 elseif isscalar(x), x = x*ones(size(n));
106 c = .01265 + .1757./n;
110 p = [255.7844, -1116.360, 1950.646, -1705.091, 745.2337, -130.2137];
111 g3 = polyval(p,x(idx));
115 idx = (x < 0.8 & x > c);
117 p = [1.91864, -8.259, 14.458, -14.6538, 6.54034, -.00022633];
120 g2 = polyval(p,(x(idx)-c1)./(.8-c1));
121 y(idx) = (.04213 + .01365*n1).*n1 .* g2;
128 g1 = sqrt(x1).*(1-x1).*(49*x1-102);
129 y(idx) = ((.0037*n1+.00078).*n1+.00006).*n1 .* g1;