C program for Kolmogorov's distribution
George Marsaglia
used the then powerful computing facilities
(SWAC) at the Institute for Numerical Analysis
to compute tables of the distribution of
Kolmogorov's goodness-of-fit test---the
maximum deviation between the true and a
sample CDF based on N sorted values from the
distribution. The calculations were based on
recursions given by Kolmogorov in his original
paper in an Italian actuarial journal in 1933.
Who knows where C or Fortran programs are available for
evaluating Kolmogorov's distribution function for given N?
(not the usual K+ or K- of Smirnov, which are relatively
easy; Kolmogorov's is max(K+,K-), or the
asymptotic form, which is also relatively easy.)
George Marsaglia
Dann Corbit
news:fWqh9.71040$5r1.2...@bin5.nnrp.aus1.giganews.com...
http://root.cern.ch/root/html/src/TMath.cxx.html
George Marsaglia
Dann Corbit <dco...@connx.com> wrote in message news:am8jq...@enews3.newsguy.com...
Unfortunately, only the asymptotic form is given there.
gm
Dann Corbit
[snip]
> > http://root.cern.ch/root/html/src/TMath.cxx.html
> >
>
> Unfortunately, only the asymptotic form is given there.
Quite right. Sorry, I did not read your original post carefully enough.
If you can point me to a PDF or PS technical paper that describes the
algorithm, I could probably write one easily enough.
--
C-FAQ: http://www.eskimo.com/~scs/C-faq/top.html
"The C-FAQ Book" ISBN 0-201-84519-9
C.A.P. FAQ: ftp://cap.connx.com/pub/Chess%20Analysis%20Project%20FAQ.htm
Dann Corbit
kl.F
A translation to C:
/*
** Crude, crappy cleanup of f2c output of CERN routine probkl.f
** D. Corbit did the cleanup. Any bugs probably belong to him.
*/
#include <math.h>
double probkl(double x)
{
static const double fj[4] = {-2., -8., -18., -32.};
int j,
i1,
i2,
i3;
double retval,
r1;
double p,
r[4],
u,
v;
u = dabs(x);
if (u < .2) {
p = 1.;
} else if (u < .755) {
/* Computing 2nd power */
r1 = u;
v = 1 / (r1 * r1);
p = 1 - (exp(v * -1.233700547316753) +
exp(v * -11.103304925850777) +
exp(v * -30.842513682918828)) * 2.50662827 / u;
} else if (u < 6.8116) {
r[1] = 0.;
r[2] = 0.;
r[3] = 0.;
/* Computing 2nd power */
r1 = u;
v = r1 * r1;
/* Computing MAX */
r1 = 3 / u;
/* I think that r1 can never be negative, so this could be simpler
*/
i2 = 1, i3 = (int) (r1 < 0 ? r1 - 0.5 : r1 + 0.5);
i1 = i2 > i3 ? i2 : i3;
for (j = 0; j < i1; ++j) {
r[j] = exp(fj[j] * v);
}
p = (r[0] - r[1] + r[2] - r[3]) * 2;
} else {
p = 0.;
}
retval = p;
return retval;