/* 
EDF-based goodness-of-fit statistics.

Last modified on 08/28/01 by Sergej V. Aksenov. 
*/

#include <math.h>

/* Macros sqr() and max() functions for double argument. */

static double sqrarg;
#define dsqr(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

static double dmaxarg1,dmaxarg2;
#define dmax(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
        (dmaxarg1) : (dmaxarg2))

/* Kolmogorov-Smirnov. */

int statistic0(double edf[],double cdf[],long n,double *result)
    {
      	int j;
      	double mdiff1, mdiff2, cdiff1, cdiff2;
 
      	mdiff1 = 0.e0;
      	mdiff2 = 0.e0;
  
      	for (j = 0; j <= n-1; j++)
	    {
		cdiff1 = edf[j] - cdf[j];
		cdiff2 = cdf[j] - edf[j] + 1.e0 / n;
		mdiff1 = dmax(cdiff1,mdiff1);
		mdiff2 = dmax(cdiff2,mdiff2);
	    }

      	*result = dmax(mdiff1,mdiff2);
      	return 0;
    }

/* Kuiper. */

int statistic1(double edf[],double cdf[],long n,double *result)
    {
      	int j;
      	double mdiff1, mdiff2, cdiff1, cdiff2;
 
      	mdiff1 = 0.e0;
      	mdiff2 = 0.e0;  

      	for (j = 0; j <= n-1; j++)
	    {
		cdiff1 = edf[j] - cdf[j];
		cdiff2 = cdf[j] - edf[j] + 1.e0 / n;
		mdiff1 = dmax(cdiff1,mdiff1);
		mdiff2 = dmax(cdiff2,mdiff2);
	    }

      	*result = mdiff1 + mdiff2;
      	return 0;
    }

/* Cramer-von Mises. */

int statistic2(double edf[],double cdf[],long n,double *result)
    {
      	int j;
      	double sum; 

      	sum = 0.e0;  
      	for (j = 0; j<= n-1; j++) sum += dsqr(cdf[j] - edf[j] +0.5 / n);

      	*result = sum + 1.0 / 12.0 / n;
      	return 0;
    }

/* Watson. */

int statistic3(double edf[],double cdf[],long n,double *result)
    {
      	int j;
      	double sum1, sum2;
 
      	sum1 = 0.e0;  
      	sum2 = 0.e0;

      	for (j = 0; j <= n-1; j++)
	    {
		sum1 += dsqr(cdf[j] - edf[j] + 0.5 / n);
		sum2 += cdf[j];
	    }  

      	*result = sum1 + 1.0 / 12.0 / n - n * dsqr(sum2 / n - 0.5);
     	return 0;
    }

/* Anderson-Darling. */

int statistic4(double edf[],double cdf[],long n,double *result)
    {
      	int j;
      	double sum; 

      	sum = 0.e0;  
      	for (j = 0; j <= n-1; j++) sum += (2*edf[j] - 1.e0/n)*(log(cdf[j]) + log(1.e0-cdf[n-1-j]));

      	*result= -n-sum;
      	return 0;
    }
	
/* Root mean squared error. */

int statistic5(double edf[],double cdf[],long n,double *result)
    {
      	int j;
      	double sum; 

      	sum = 0.e0;  
      	for (j = 0; j <= n-1; j++) sum += dsqr(cdf[j] - edf[j]);

      	*result=sqrt(sum);
      	return 0;
    }

#undef dsqr
#undef dmax