Evaluating the Anderson-darling Distribution

Except for n = 1, only the limit as n → ∞ for the distribution of the Anderson-Darling test for uniformity has been found, and that in so complicated a form that published values for a few percentiles had to be determined by numerical integration, saddlepoint or other approximation methods. We give here our method for evaluating that asymptotic distribution to great accuracy—directly, via series with two-term recursions. We also give, for any particular n, a procedure for evaluating the distribution to the fourth digit, based on empirical CDF's from samples of size 10 10 .

This is a special case of the Cramer-von Mises approach: use 1 0 w(x)[Fn(x)−x] 2 dx, the squared area between the sample CDF Fn(x) (the staircase) and the diagonal y = x, using a suitable weight function w(x).That approach: use the weighted square of the area, is in contrast to the Kolmogorov approach: use the maximum distance from the staircase to the diagonal, (for which the distribution has only recently become available [6]).Choice of the weight function w(x) = 1 x(1−x) was suggested by L.J. Savage, as it divides the distances from the diagonal to the corners of the staircase by their standard deviations.That choice leads to the statistic An above, for which Anderson and Darling [1] derived a complicated expression for the asymptotic distribution.
Thus, with Fn(x) the sample CDF of a set {X1, ..., Xn} of iid random variables in [0,1), Fn(x) = number of X1, X2, . . ., Xn that are ≤ x n , the Anderson-Darling statistic for testing that the X's came from a uniform distribution is Since Fn(x) is a step function, expressing the above integral as a sum in two parts leads to a collection of elementary integrals from which a little manipulation provides the form Even for n = 2 the distribution is difficult to evaluate (via numerical integration), and for specific n > 2, all that seems available is some tabled values for n ≤ 8 based on simulations by Lewis [4] -simulations limited by CPU speeds and sampling methods circa 1960.But current speeds of CPU's, and fast methods for generating an ordered sample of uniform variates by means of the ziggurat method for exponential variates, make it feasible to generate 10 10 random values of An for n's in the hundreds, and thus determine the distribution with considerable accuracy.We do this for n = 8, 16, 32, 64, 128 and the results lead to formulas for in-between n's that seem to provide Pr(An < z) with accuracies to the fifth digit, inferred from the size of the samples used to derive them and supported by extensive testing.The limiting distribution of An also suffers from lack of a method for its accurate determination.Anderson,Darling [2] and Lewis [4] both give the same values for the 90,95 and 99 percentiles (their 99th percentile, 3.857 is actually 3.878125...).Other approximations are given by Sinclair and Spurr [7], while Giles [3] has recently suggested a saddlepoint method.
We provide here a method for determining that asymptotic distribution to arbitrary precision, with a C version to double precision accuracy, roughly 13-15 digits.
We first describe the method for evaluating the asymptotic distribution, limn→∞ Pr(An < z), to desired accuracy, and then show how it may be used to determine, for given n, Pr(An < z) by means of adjustments based on simulations of size 10 10 .
2 The limiting distribution of A n .
An expression for the limiting distribution of An was given by Anderson and Darling [1].The method was based on a development of Doob for the absorption probability of a diffusion model.They gave This is a strange distribution function.Anderson and Darling [2] used numerical integration to find the 90, 95 and 99 percentiles.(They are reported as 1.933,2.492and 3.857; the true values to 20 places are 1.9329578327415937304, 2.4923671600494096176 and 3.8781250216053948842.) Lewis [4], also using numerical integration, published a table giving lim Pr(An < z) with 4-place accuracy for selected z values, as well as the same three percentiles with the wrong 3.857 value for 3.878125.... Other values have been provided by Sinclair and Spurr [6], (approximate inversion of the characterisitic function), and Giles [3], (saddlepoint approximations).In all, it seems that relatively few values or percentiles have been provided, all by approximation methods and sometimes giving less than the claimed 3-4 digits of accuracy.Note that Sinclair and Spur report a better value, 3.880 as the 99 percentile 3.878125..., which Giles disputes in a footnote, sticking to 3.857 as the 'true' value, presumably because it was given by both Anderson-Darling [2] and Lewis [4].
We will provide a method for evaluating the above distribution with accuracy limited to the computer's ability to distinguish between floating point numbers, give a C program for implementing it, and also give a quick-and-easy approximation that gives accuracy better than .000002for probabilities less than .9and .0000008for those beyond.
We designate the limiting distribution by ADinf(z), and put it in the form where The difficulty is in evaluating f (z, j).To do that, we expand exp( ) in a series: with, for given j, and tj, We then have, all with fixed j and t = tj = (4j + 1) 2 π 2 /(8z): and-the key to accurate evaluation of ADinf()-the recursion: Using these recursions, for fixed j in the series for f (z, j), we may then evaluate ADinf(z) as a series, The accompanying C version will evaluate ADinf(z) to around fifteen places, checked with a 30-digit Maple version.Note that for the two initial values of the recursions, the second requires erfc, the complementary error function.Our version of the complementary normal distribution function, cPhi(x), is included and recommended for use, as it is more accurate than many of the available erf(x) or erfc(x) routines in C compiler libraries.
Figure 2 shows the distribution and density of A∞.The distribution is infinitely flat (every derivative is zero) at z = 0 and nearly flat as z passes 6 or so, with a slow approach to 1.The C (or Maple) procedure gives ADinf(9) = .999960465988611(.999960465988612484992562014458), and ADinf(10) = .999986184964588(.999986184964589314168018038088).The above method for evaluating ADinf(z) to unlimited accuracy, requiring the normal integral and a loop within a loop, might be more complicated and/or have more power than many users want to invokecertainly in applications where 6-digit accuracy may be more than adequate.For that reason we offer the following two-piece formula, with accuracies indicated, (using lowercase adinf rather than the mixed ADinf): The function ADinf(z) starts like 2z −1/2 e π 2 /(8z) , and π 2 /8 = 1.23370055 is changed to 1.2337141 to ensure (7-place) continutity at z = 2, with a fifth degree polynomial in Horner form as a multiplier for the range 0 < z ≤ 2. With increasing z, we might expect a standard extreme-value form: e −e a−bz , and examination shows that ln(− ln(ADinf(z))) looks quite linear.But a linear form in that top exponent does not provide the accuracy we seek, and a Horner polynomial to degree 5 is used.
3 Getting the distribution of A n by simulation.
Until now, little was known about the distribution of An for various n.Lewis [4] gave tables for small n, based on simulations conducted at the IBM Research Center, but was unable to go beyond n = 8.He concluded that the convergence of the distribution of An to that of A∞ is "quite rapid".Such confidence in a rapid approach seems to hold only for probabilities greater than about .8.We think that use of A∞'s distribution for that of An is not good enough for modern computer-aided statistical procedures.The worst error in using A∞ for a particular An is about .044/n,near the 33rd percentile.But there is some support for Lewis's confidence: the approach is more rapid than that for the limiting form of Kolmogorov's distribution, which has worst error of about .278/√ n, also around the 33rd percentile [6].
By developing very fast ways to generate a sequence of ordered uniform variables, we will be able to get samples of An large enough (typically 10 10 with each of n = 8, 16, 32, 64, 128) to provide accuracies better than .00005for determining Pr(An < z) or for converting an observed An to a p-value.
To do this, we need a way to partition the possible values of An into cells with approximately equal probabilities.But this need not be done directly.When applying a goodness-of-fit test to the values x1 < x2 < • • • < xn, whatever version of the Kolmogorov-Smirnov class of tests we use, it is necessary to convert the resulting statistic to a uniform variable, a p-value, in [0,1).We do not need Pr(An < z) directly, but, rather, we need a way to convert the observed value to a uniform [0,1) value.
We do that in the following way: Since the distribution of An is close to that of A∞, if Z is a random value with distribution of An, then ADinf(Z) should be close to uniformly distributed in [0,1).Therefore, if we divide the unit interval into 1000 cells, 0 to .001,.001 to .002, . ..,.999 to 1, then count the number of times that ADinf(Z) falls into those intervals.We need only make a small adjustment to make the resulting empirical distribution quite close to uniform in [0,1)-provided our mean cell counts are large enough to provide the necessary accuracy.Our mean cell counts are around 10 7 .
A computer approximation to the true errfix(n, x) needs to have sufficient accuracy to ensure that the composite result, ADinf(z) + errfix(n, ADinf(z)) will be close enough to uniform for practical applications.
The error function errfix(n, x), 0 < x < 1 given below is based on extensive simulation.For a given value of the random variable Z = An, we evaluate x = ADinf(Z), then convert that x to a uniform [0,1) random variable p by means of p = x + errfix(n, x).Graphs of errfix(n, x), 0 < x < 1 are given in Figure 3, for n = 8, 16, 32, 64, 128.Here, we get lucky.We are able to adequately represent, piecewise, any one of those error curves by means of three fixed functions g1, g2, g3, adjusting the scale as a function of 1/n and providing stretch factors in the g's arguments.Each error curve is made of three parts, an initial dip that reaches the x-axis at a point c = c(n) and can be adequately represented as g1(x/c) times a function of 1/n.Then the second region from c to .8 can be adequately represented as a function of 1/n times g2(x − c)/(.8 − c), and finally, a last dip from .8 to 1 than can be represented as g3(x)/n.For each error curve, the right crossing is satisfactorily close to x = .8for all n, but the left crossing, c, varies enough to require an empirical function of n: c(n) = .01265+ .1757/n.
Since the exact distribution of Z = An is not known, but only approximated from samples of 10 10 for various n, we base our claim on:

Figure 1 :
Figure 1: The density and distribution function of A ∞ .