summaryrefslogtreecommitdiff
path: root/libm/double/kolmogorov.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/kolmogorov.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/double/kolmogorov.c')
-rw-r--r--libm/double/kolmogorov.c243
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