lib/libm/kolmogorov.c

   1
   2 /* Re Kolmogorov statistics, here is Birnbaum and Tingey's formula for the
   3    distribution of D+, the maximum of all positive deviations between a
   4    theoretical distribution function P(x) and an empirical one Sn(x)
   5    from n samples.
   6
   7      +
   8     D  =         sup     [P(x) - S (x)]
   9      n     -inf < x < inf         n
  10
  11
  12                   [n(1-e)]
  13         +            -                    v-1              n-v
  14     Pr{D   > e} =    >    C    e (e + v/n)    (1 - e - v/n)
  15         n            -   n v
  16                     v=0
  17
  18     [n(1-e)] is the largest integer not exceeding n(1-e).
  19     nCv is the number of combinations of n things taken v at a time.  */
  20
  21
  22 #include "mconf.h"
  23 #ifdef ANSIPROT
  24 extern double pow ( double, double );
  25 extern double floor ( double );
  26 extern double lgam ( double );
  27 extern double exp ( double );
  28 extern double sqrt ( double );
  29 extern double log ( double );
  30 extern double fabs ( double );
  31 double smirnov ( int, double );
  32 double kolmogorov ( double );
  33 #else
  34 double pow (), floor (), lgam (), exp (), sqrt (), log (), fabs ();
  35 double smirnov (), kolmogorov ();
  36 #endif
  37 extern double MAXLOG;
  38
  39 /* Exact Smirnov statistic, for one-sided test.  */
  40 double
  41 smirnov (n, e)
  42      int n;
  43      double e;
  44 {
  45   int v, nn;
  46   double evn, omevn, p, t, c, lgamnp1;
  47
  48   if (n <= 0 || e < 0.0 || e > 1.0)
  49     return (-1.0);
  50   nn = floor ((double) n * (1.0 - e));
  51   p = 0.0;
  52   if (n < 1013)
  53     {
  54       c = 1.0;
  55       for (v = 0; v <= nn; v++)
  56         {
  57           evn = e + ((double) v) / n;
  58           p += c * pow (evn, (double) (v - 1))
  59             * pow (1.0 - evn, (double) (n - v));
  60           /* Next combinatorial term; worst case error = 4e-15.  */
  61           c *= ((double) (n - v)) / (v + 1);
  62         }
  63     }
  64   else
  65     {
  66       lgamnp1 = lgam ((double) (n + 1));
  67       for (v = 0; v <= nn; v++)
  68         {
  69           evn = e + ((double) v) / n;
  70           omevn = 1.0 - evn;
  71           if (fabs (omevn) > 0.0)
  72             {
  73               t = lgamnp1
  74                 - lgam ((double) (v + 1))
  75                 - lgam ((double) (n - v + 1))
  76                 + (v - 1) * log (evn)
  77                 + (n - v) * log (omevn);
  78               if (t > -MAXLOG)
  79                 p += exp (t);
  80             }
  81         }
  82     }
  83   return (p * e);
  84 }
  85
  86
  87 /* Kolmogorov's limiting distribution of two-sided test, returns
  88    probability that sqrt(n) * max deviation > y,
  89    or that max deviation > y/sqrt(n).
  90    The approximation is useful for the tail of the distribution
  91    when n is large.  */
  92 double
  93 kolmogorov (y)
  94      double y;
  95 {
  96   double p, t, r, sign, x;
  97
  98   x = -2.0 * y * y;
  99   sign = 1.0;
 100   p = 0.0;
 101   r = 1.0;
 102   do
 103     {
 104       t = exp (x * r * r);
 105       p += sign * t;
 106       if (t == 0.0)
 107         break;
 108       r += 1.0;
 109       sign = -sign;
 110     }
 111   while ((t / p) > 1.1e-16);
 112   return (p + p);
 113 }
 114
 115 /* Functional inverse of Smirnov distribution
 116    finds e such that smirnov(n,e) = p.  */
 117 double
 118 smirnovi (n, p)
 119      int n;
 120      double p;
 121 {
 122   double e, t, dpde;
 123
 124   if (p <= 0.0 || p > 1.0)
 125     {
 126       mtherr ("smirnovi", DOMAIN);
 127       return 0.0;
 128     }
 129   /* Start with approximation p = exp(-2 n e^2).  */
 130   e = sqrt (-log (p) / (2.0 * n));
 131   do
 132     {
 133       /* Use approximate derivative in Newton iteration. */
 134       t = -2.0 * n * e;
 135       dpde = 2.0 * t * exp (t * e);
 136       if (fabs (dpde) > 0.0)
 137         t = (p - smirnov (n, e)) / dpde;
 138       else
 139         {
 140           mtherr ("smirnovi", UNDERFLOW);
 141           return 0.0;
 142         }
 143       e = e + t;
 144       if (e >= 1.0 || e <= 0.0)
 145         {
 146           mtherr ("smirnovi", OVERFLOW);
 147           return 0.0;
 148         }
 149     }
 150   while (fabs (t / e) > 1e-10);
 151   return (e);
 152 }
 153
 154
 155 /* Functional inverse of Kolmogorov statistic for two-sided test.
 156    Finds y such that kolmogorov(y) = p.
 157    If e = smirnovi (n,p), then kolmogi(2 * p) / sqrt(n) should
 158    be close to e.  */
 159 double
 160 kolmogi (p)
 161      double p;
 162 {
 163   double y, t, dpdy;
 164
 165   if (p <= 0.0 || p > 1.0)
 166     {
 167       mtherr ("kolmogi", DOMAIN);
 168       return 0.0;
 169     }
 170   /* Start with approximation p = 2 exp(-2 y^2).  */
 171   y = sqrt (-0.5 * log (0.5 * p));
 172   do
 173     {
 174       /* Use approximate derivative in Newton iteration. */
 175       t = -2.0 * y;
 176       dpdy = 4.0 * t * exp (t * y);
 177       if (fabs (dpdy) > 0.0)
 178         t = (p - kolmogorov (y)) / dpdy;
 179       else
 180         {
 181           mtherr ("kolmogi", UNDERFLOW);
 182           return 0.0;
 183         }
 184       y = y + t;
 185     }
 186   while (fabs (t / y) > 1e-10);
 187   return (y);
 188 }
 189
 190
 191 #ifdef SALONE
 192 /* Type in a number.  */
 193 void
 194 getnum (s, px)
 195      char *s;
 196      double *px;
 197 {
 198   char str[30];
 199
 200   printf (" %s (%.15e) ? ", s, *px);
 201   gets (str);
 202   if (str[0] == '\0' || str[0] == '\n')
 203     return;
 204   sscanf (str, "%lf", px);
 205   printf ("%.15e\n", *px);
 206 }
 207
 208 /* Type in values, get answers.  */
 209 void
 210 main ()
 211 {
 212   int n;
 213   double e, p, ps, pk, ek, y;
 214
 215   n = 5;
 216   e = 0.0;
 217   p = 0.1;
 218 loop:
 219   ps = n;
 220   getnum ("n", &ps);
 221   n = ps;
 222   if (n <= 0)
 223     {
 224       printf ("? Operator error.\n");
 225       goto loop;
 226     }
 227   /*
 228   getnum ("e", &e);
 229   ps = smirnov (n, e);
 230   y = sqrt ((double) n) * e;
 231   printf ("y = %.4e\n", y);
 232   pk = kolmogorov (y);
 233   printf ("Smirnov = %.15e, Kolmogorov/2 = %.15e\n", ps, pk / 2.0);
 234 */
 235   getnum ("p", &p);
 236   e = smirnovi (n, p);
 237   printf ("Smirnov e = %.15e\n", e);
 238   y = kolmogi (2.0 * p);
 239   ek = y / sqrt ((double) n);
 240   printf ("Kolmogorov e = %.15e\n", ek);
 241   goto loop;
 242 }
 243 #endif