diff options
Diffstat (limited to 'libm/double/kolmogorov.c')
| -rw-r--r-- | libm/double/kolmogorov.c | 243 |
1 files changed, 0 insertions, 243 deletions
diff --git a/libm/double/kolmogorov.c b/libm/double/kolmogorov.c deleted file mode 100644 index 0d6fe92bd..000000000 --- a/libm/double/kolmogorov.c +++ /dev/null @@ -1,243 +0,0 @@ - -/* 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 <math.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 |