computed Kolmogorov-Smirnov's D

diff --git a/kmath.c b/kmath.c
index b6fbf60f97012dce8af5f8ec10373f93ff5393bc..9807b00ee5fee69a3903c734b87bddd2af9039f8 100644 (file)
--- a/kmath.c
+++ b/kmath.c
@@ -416,12 +416,37 @@ double kf_betai(double a, double b, double x)
        return x < (a + 1.) / (a + b + 2.)? kf_betai_aux(a, b, x) : 1. - kf_betai_aux(b, a, 1. - x);
 }
 
+/******************
+ *** Statistics ***
+ ******************/
+
+double km_ks_dist(int na, const double a[], int nb, const double b[]) // a[] and b[] MUST BE sorted
+{
+       int ia = 0, ib = 0;
+       double fa = 0, fb = 0, sup = 0, na1 = 1. / na, nb1 = 1. / nb;
+       while (ia < na || ib < nb) {
+               if (ia == na) fb += nb1, ++ib;
+               else if (ib == nb) fa += na1, ++ia;
+               else if (a[ia] < b[ib]) fa += na1, ++ia;
+               else if (a[ia] > b[ib]) fb += nb1, ++ib;
+               else fa += na1, fb += nb1, ++ia, ++ib;
+               if (sup < fabs(fa - fb)) sup = fabs(fa - fb);
+       }
+       return sup;
+}
+
 #ifdef KF_MAIN
 #include <stdio.h>
+#include "ksort.h"
+KSORT_INIT_GENERIC(double)
 int main(int argc, char *argv[])
 {
        double x = 5.5, y = 3;
        double a, b;
+       double xx[] = {0.22, -0.87, -2.39, -1.79, 0.37, -1.54, 1.28, -0.31, -0.74, 1.72, 0.38, -0.17, -0.62, -1.10, 0.30, 0.15, 2.30, 0.19, -0.50, -0.09};
+       double yy[] = {-5.13, -2.19, -2.43, -3.83, 0.50, -3.25, 4.32, 1.63, 5.18, -0.43, 7.11, 4.87, -3.10, -5.81, 3.76, 6.31, 2.58, 0.07, 5.76, 3.50};
+       ks_introsort(double, 20, xx); ks_introsort(double, 20, yy);
+       printf("K-S distance: %f\n", km_ks_dist(20, xx, 20, yy));
        printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
        printf("upper-gamma(%lg,%lg): %lg\n", x, y, kf_gammaq(y, x)*tgamma(y));
        a = 2; b = 2; x = 0.5;