https://github.com/ntamas/plfit
kolmogorov.c
/* kolmogorov.c
*
* Copyright (C) 2010-2011 Tamas Nepusz
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <math.h>
#include "kolmogorov.h"
double plfit_kolmogorov(double z) {
const double fj[4] = { -2, -8, -18, -32 };
const double w = 2.50662827;
const double c1 = -1.2337005501361697; /* -pi^2 / 8 */
const double c2 = -11.103304951225528; /* 9*c1 */
const double c3 = -30.842513753404244; /* 25*c1 */
double u = fabs(z);
double v;
if (u < 0.2)
return 1;
if (u < 0.755) {
v = 1.0 / (u*u);
return 1 - w * (exp(c1*v) + exp(c2*v) + exp(c3*v)) / u;
}
if (u < 6.8116) {
double r[4] = { 0, 0, 0, 0 };
long int maxj = (long int)(3.0 / u + 0.5);
long int j;
if (maxj < 1)
maxj = 1;
v = u*u;
for (j = 0; j < maxj; j++) {
r[j] = exp(fj[j] * v);
}
return 2*(r[0] - r[1] + r[2] - r[3]);
}
return 0;
}
double plfit_ks_test_one_sample_p(double d, size_t n) {
return plfit_kolmogorov(d * sqrt((double) n));
}
double plfit_ks_test_two_sample_p(double d, size_t n1, size_t n2) {
return plfit_kolmogorov(d * sqrt(n1*n2 / ((double)(n1+n2))));
}