summaryrefslogtreecommitdiff
path: root/libm/double/kolmogorov.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/kolmogorov.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/kolmogorov.c')
-rw-r--r--libm/double/kolmogorov.c243
1 files changed, 243 insertions, 0 deletions
diff --git a/libm/double/kolmogorov.c b/libm/double/kolmogorov.c
new file mode 100644
index 000000000..0d6fe92bd
--- /dev/null
+++ b/libm/double/kolmogorov.c
@@ -0,0 +1,243 @@
+
+/* 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