Quelle  kolmogorov.c

/* Re Kolmogorov statistics, here is Birnbaum and Tingey's formula for the
   distribution of D+, the maximum of all positive deviations between a
   theoretical distribution function P(x) and an empirical one Sn(x)
   from n samples.

     +
    D  =         sup     [P(x) - S (x)]
     n     -inf < x < inf         n


                  [n(1-e)]
        +            -                    v-1              n-v
    Pr{D   > e} =    >    C    e (e + v/n)    (1 - e - v/n)
        n            -   n v
                    v=0

    [n(1-e)] is the largest integer not exceeding n(1-e).
    nCv is the number of combinations of n things taken v at a time.  */

#include "mconf.h"
#ifdef ANSIPROT
extern double pow(double, double);
extern double floor(double);
extern double lgam(double);
extern double exp(double);
extern double sqrt(double);
extern double log(double);
extern double fabs(double);
double smirnov(int, double);
double kolmogorov(double);
#else
double pow(), floor(), lgam(), exp(), sqrt(), log(), fabs();
double smirnov(), kolmogorov();
#endif
extern double MAXLOG;

/* Exact Smirnov statistic, for one-sided test.  */
double smirnov(n, e) int n;
double e;
{
  int v, nn;
  double evn, omevn, p, t, c, lgamnp1;

  if (n <= 0 || e < 0.0 || e > 1.0)
    return (-1.0);
  nn = floor((double)n * (1.0 - e));
  p = 0.0;
  if (n < 1013) {
    c = 1.0;
    for (v = 0; v <= nn; v++) {
      evn = e + ((double)v) / n;
      p += c * pow(evn, (double)(v - 1)) * pow(1.0 - evn, (double)(n - v));
      /* Next combinatorial term; worst case error = 4e-15.  */
      c *= ((double)(n - v)) / (v + 1);
    }
  } else {
    lgamnp1 = lgam((double)(n + 1));
    for (v = 0; v <= nn; v++) {
      evn = e + ((double)v) / n;
      omevn = 1.0 - evn;
      if (fabs(omevn) > 0.0) {
        t = lgamnp1 - lgam((double)(v + 1)) - lgam((double)(n - v + 1)) +
            (v - 1) * log(evn) + (n - v) * log(omevn);
        if (t > -MAXLOG)
          p += exp(t);
      }
    }
  }
  return (p * e);
}

/* Kolmogorov's limiting distribution of two-sided test, returns
   probability that sqrt(n) * max deviation > y,
   or that max deviation > y/sqrt(n).
   The approximation is useful for the tail of the distribution
   when n is large.  */
double kolmogorov(y) double y;
{
  double p, t, r, sign, x;

  x = -2.0 * y * y;
  sign = 1.0;
  p = 0.0;
  r = 1.0;
  do {
    t = exp(x * r * r);
    p += sign * t;
    if (t == 0.0)
      break;
    r += 1.0;
    sign = -sign;
  } while ((t / p) > 1.1e-16);
  return (p + p);
}

/* Functional inverse of Smirnov distribution
   finds e such that smirnov(n,e) = p.  */
double smirnovi(n, p) int n;
double p;
{
  double e, t, dpde;

  if (p <= 0.0 || p > 1.0) {
    mtherr("smirnovi", DOMAIN);
    return 0.0;
  }
  /* Start with approximation p = exp(-2 n e^2).  */
  e = sqrt(-log(p) / (2.0 * n));
  do {
    /* Use approximate derivative in Newton iteration. */
    t = -2.0 * n * e;
    dpde = 2.0 * t * exp(t * e);
    if (fabs(dpde) > 0.0)
      t = (p - smirnov(n, e)) / dpde;
    else {
      mtherr("smirnovi", UNDERFLOW);
      return 0.0;
    }
    e = e + t;
    if (e >= 1.0 || e <= 0.0) {
      mtherr("smirnovi", OVERFLOW);
      return 0.0;
    }
  } while (fabs(t / e) > 1e-10);
  return (e);
}

/* Functional inverse of Kolmogorov statistic for two-sided test.
   Finds y such that kolmogorov(y) = p.
   If e = smirnovi (n,p), then kolmogi(2 * p) / sqrt(n) should
   be close to e.  */
double kolmogi(p) double p;
{
  double y, t, dpdy;

  if (p <= 0.0 || p > 1.0) {
    mtherr("kolmogi", DOMAIN);
    return 0.0;
  }
  /* Start with approximation p = 2 exp(-2 y^2).  */
  y = sqrt(-0.5 * log(0.5 * p));
  do {
    /* Use approximate derivative in Newton iteration. */
    t = -2.0 * y;
    dpdy = 4.0 * t * exp(t * y);
    if (fabs(dpdy) > 0.0)
      t = (p - kolmogorov(y)) / dpdy;
    else {
      mtherr("kolmogi", UNDERFLOW);
      return 0.0;
    }
    y = y + t;
  } while (fabs(t / y) > 1e-10);
  return (y);
}

#ifdef SALONE
/* Type in a number.  */
void getnum(s, px) char *s;
double *px;
{
  char str[30];

  printf(" %s (%.15e) ? ", s, *px);
  gets(str);
  if (str[0] == '\0' || str[0] == '\n')
    return;
  sscanf(str, "%lf", px);
  printf("%.15e\n", *px);
}

/* Type in values, get answers.  */
void main() {
  int n;
  double e, p, ps, pk, ek, y;

  n = 5;
  e = 0.0;
  p = 0.1;
loop:
  ps = n;
  getnum("n", &ps);
  n = ps;
  if (n <= 0) {
    printf("? Operator error.\n");
    goto loop;
  }
  /*
  getnum ("e", &e);
  ps = smirnov (n, e);
  y = sqrt ((double) n) * e;
  printf ("y = %.4e\n", y);
  pk = kolmogorov (y);
  printf ("Smirnov = %.15e, Kolmogorov/2 = %.15e\n", ps, pk / 2.0);
*/
  getnum("p", &p);
  e = smirnovi(n, p);
  printf("Smirnov e = %.15e\n", e);
  y = kolmogi(2.0 * p);
  ek = y / sqrt((double)n);
  printf("Kolmogorov e = %.15e\n", ek);
  goto loop;
}
#endif