sddsdistest.c File Reference

Detailed Description

Statistical distribution testing tool for SDDS datasets.

This program performs statistical tests on columns of data from an SDDS file against specified distributions (Gaussian, Poisson, Student's t, or Chi-Squared). It supports both Kolmogorov-Smirnov (KS) and Chi-Squared tests, along with the ability to use user-defined distributions provided via external files.

Usage

sddsdistest [<inputfile>] [<outputfile>]
            [-pipe=[in][,out]]
             -column=<name>[,sigma=<name>]...
            [-exclude=<name>[,...]]
            [-degreesOfFreedom={<value>|@<parameterName>}]
            [-test={ks|chisquared}]
            [-fileDistribution=<filename>,<indepName>,<depenName>]
            [-gaussian]
            [-poisson]
            [-student]
            [-chisquared]
            [-majorOrder=row|column]
            [-threads=<number>]

Options

Required	Description
-column	Specifies the columns to test, with optional sigma columns for error handling.

Optional	Description
-pipe	Enables piping for input and/or output data streams.
-exclude	Excludes specified columns from testing.
-degreesOfFreedom	Specifies degrees of freedom as a fixed value or a parameter.
-test	Selects the test to perform: Kolmogorov-Smirnov or Chi-Squared.
-fileDistribution	Uses a user-defined distribution from a file.
-gaussian	Tests against a Gaussian distribution.
-poisson	Tests against a Poisson distribution.
-student	Tests against a Student's t distribution.
-chisquared	Tests against a Chi-Squared distribution.
-majorOrder	Specifies the data ordering: 'row' or 'column'.
-threads	Number of threads for parallel computations.

Incompatibilities

• The -test option is incompatible with specifying multiple distributions (e.g., -gaussian and -poisson).
• Only one of the following may be specified:
  • -fileDistribution
  • -gaussian
  • -poisson
  • -student
  • -chisquared

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, R. Soliday, Xuesong Jiao, H. Shang

Definition in file sddsdistest.c.

#include "mdb.h"
#include "SDDS.h"
#include "SDDSutils.h"
#include "scan.h"

Go to the source code of this file.

Functions
void	compareToFileDistribution (char *output, long testCode, SDDS_DATASET *SDDSin, char **columnName, long columnNames, char *distFile, char *distFileIndep, char *distFileDepen)
void	compareToDistribution (char *output, long testCode, SDDS_DATASET *SDDSin, char **columnName, char **sigmaName, long columnNames, long distCode, long degreesFree, char *dofParameter, short columnMajorOrder, int threads)
void	ksTestWithFunction (double *data, int64_t rows, double(*CDF)(double x), double *statReturn, double *sigLevelReturn)
void	chiSquaredTestWithFunction (double *data, int64_t rows, double(*PDF)(double x), double *statReturn, double *sigLevelReturn)
int	main (int argc, char **argv)
double	gaussianPDF (double x)
double	gaussianCDF (double x)
double	poissonPDF (double xd)
double	poissonCDF (double xd)
double	studentPDF (double t)
double	studentCDF (double t)
double	chiSquaredPDF (double x)
double	chiSquaredCDF (double x)

Function Documentation

chiSquaredCDF()

double chiSquaredCDF

(

double

x

)

Definition at line 402 of file sddsdistest.c.

                               {
  double chiSqr;
 
  if (x < 0)
    x = 0;
  chiSqr = x * DOF / sampleMean;
  return 1 - gammaQ(DOF / 2.0, chiSqr / 2.0);
}

chiSquaredPDF()

double chiSquaredPDF

(

double

x

)

Definition at line 392 of file sddsdistest.c.

                               {
  double chiSqr, DOFover2;
 
  if (x < 0)
    return 0;
  chiSqr = x * DOF / sampleMean;
  DOFover2 = DOF / 2.0;
  return exp((DOFover2 - 1.0) * log(chiSqr) - chiSqr / 2 - DOFover2 * LOG2 - lgamma(DOFover2));
}

chiSquaredTestWithFunction()

void chiSquaredTestWithFunction	(	double *	data,
		int64_t	rows,
		double(*	PDF )(double x),
		double *	statReturn,
		double *	sigLevelReturn )

Definition at line 532 of file sddsdistest.c.

                                                                                                                                 {
  SDDS_Bomb("Chi-squared distribution test not implemented yet---wouldn't you really like a nice K-S test instead?");
}

compareToDistribution()

void compareToDistribution	(	char *	output,
		long	testCode,
		SDDS_DATASET *	SDDSin,
		char **	columnName,
		char **	sigmaName,
		long	columnNames,
		long	distCode,
		long	degreesFree,
		char *	dofParameter,
		short	columnMajorOrder,
		int	threads )

Definition at line 411 of file sddsdistest.c.

                                                                                                                                                                                                                               {
  SDDS_DATASET SDDSout;
  double *data, *data1, stat, sigLevel;
  long iStat, iSigLevel, iCount, iColumnName;
  long icol;
  int64_t rows, row;
  iStat = iSigLevel = iCount = iColumnName = 0;
  if (!SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, "sddsdistest output", output) ||
      0 > SDDS_DefineParameter(&SDDSout, "distestDistribution", NULL, NULL, "sddsdistest distribution name", NULL,
                               SDDS_STRING, option[distCode]) ||
      0 > SDDS_DefineParameter(&SDDSout, "distestTest", NULL, NULL, "sddsdistest test name", NULL, SDDS_STRING, testChoice[testCode]) ||
      0 > (iCount = SDDS_DefineParameter(&SDDSout, "Count", NULL, NULL, "Number of data points", NULL, SDDS_LONG, 0))) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  if (columnMajorOrder != -1)
    SDDSout.layout.data_mode.column_major = columnMajorOrder;
  else
    SDDSout.layout.data_mode.column_major = SDDSin->layout.data_mode.column_major;
  switch (testCode) {
  case KS_TEST:
    if (0 > (iColumnName = SDDS_DefineColumn(&SDDSout, "ColumnName", NULL, NULL, "Column analysed by sddsdistest",
                                             NULL, SDDS_STRING, 0)) ||
        0 > (iStat = SDDS_DefineColumn(&SDDSout, "D", NULL, NULL, "Kolmogorov-Smirnov D statistic", NULL, SDDS_DOUBLE, 0)) ||
        0 > (iSigLevel = SDDS_DefineColumn(&SDDSout, "distestSigLevel", "P(D$ba$n>D)", NULL, "Probability of exceeding D", NULL, SDDS_DOUBLE, 0))) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    break;
  case CHI_TEST:
    if (0 > (iColumnName = SDDS_DefineColumn(&SDDSout, "ColumnName", NULL, NULL, "Column analysed by sddsdistest",
                                             NULL, SDDS_STRING, 0)) ||
        0 > (iStat = SDDS_DefineParameter(&SDDSout, "ChiSquared", "$gh$r$a2$n", NULL, "Chi-squared statistic", NULL,
                                          SDDS_DOUBLE, 0)) ||
        0 > (iSigLevel = SDDS_DefineParameter(&SDDSout, "distestSigLevel", "P($gh$r$a2$n$ba$n>$gh$r$a2$n)", NULL, "Probability of exceeding $gh$r$a2$n", NULL, SDDS_DOUBLE, 0))) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    break;
  default:
    SDDS_Bomb("Invalid testCode seen in compareToDistribution--this shouldn't happen.");
    break;
  }
  if (!SDDS_WriteLayout(&SDDSout))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  while (SDDS_ReadPage(SDDSin) > 0) {
    if (!SDDS_StartPage(&SDDSout, columnNames))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    rows = SDDS_CountRowsOfInterest(SDDSin);
    for (icol = 0; icol < columnNames; icol++) {
      stat = 0;
      sigLevel = 1;
      if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, columnName, columnNames, iColumnName) ||
          !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCount, rows, -1)) {
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
      }
      if (rows >= 2) {
        if (!(data = SDDS_GetColumnInDoubles(SDDSin, columnName[icol])))
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
        if (dofParameter)
          SDDS_GetParameterAsLong(SDDSin, dofParameter, &DOF);
        else
          DOF = degreesFree;
        switch (distCode) {
        case CLO_GAUSSIAN:
          computeMomentsThreaded(&sampleMean, NULL, &sampleStDev, NULL, data, rows, threads);
          if (testCode == KS_TEST)
            ksTestWithFunction(data, rows, gaussianCDF, &stat, &sigLevel);
          else
            chiSquaredTestWithFunction(data, rows, gaussianPDF, &stat, &sigLevel);
          break;
        case CLO_POISSON:
          computeMomentsThreaded(&sampleMean, NULL, NULL, NULL, data, rows, threads);
          if (testCode == KS_TEST)
            ksTestWithFunction(data, rows, poissonCDF, &stat, &sigLevel);
          else
            chiSquaredTestWithFunction(data, rows, poissonPDF, &stat, &sigLevel);
          break;
        case CLO_STUDENT:
          if (DOF < 1)
            SDDS_Bomb("must have at least one degree of freedom for Student distribution tests");
          computeMomentsThreaded(&sampleMean, NULL, NULL, NULL, data, rows, threads);
          if (sigmaName && sigmaName[icol]) {
            if (!(data1 = SDDS_GetColumnInDoubles(SDDSin, sigmaName[icol])))
              SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
            for (row = 0; row < rows; row++)
              /* compute t = (x-mu)/sigma */
              data[row] = (data[row] - sampleMean) / data1[row];
            free(data1);
          } else {
            for (row = 0; row < rows; row++)
              data[row] -= sampleMean;
          }
          if (testCode == KS_TEST)
            ksTestWithFunction(data, rows, studentCDF, &stat, &sigLevel);
          else
            chiSquaredTestWithFunction(data, rows, studentPDF, &stat, &sigLevel);
          break;
        case CLO_CHISQUARED:
          computeMomentsThreaded(&sampleMean, NULL, NULL, NULL, data, rows, threads);
          if (DOF < 1)
            SDDS_Bomb("must have at least one degree of freedom for chi-squared distribution tests");
          if (testCode == KS_TEST)
            ksTestWithFunction(data, rows, chiSquaredCDF, &stat, &sigLevel);
          else
            chiSquaredTestWithFunction(data, rows, chiSquaredPDF, &stat, &sigLevel);
          break;
        default:
          SDDS_Bomb("Invalid distCode in compareToDistribution--this shouldn't happen");
          break;
        }
      }
      if (!SDDS_SetRowValues(&SDDSout, SDDS_PASS_BY_VALUE | SDDS_SET_BY_INDEX, icol, iStat, stat, iSigLevel, sigLevel, -1))
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (!SDDS_WritePage(&SDDSout))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  if (!SDDS_Terminate(&SDDSout)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    exit(EXIT_FAILURE);
  }
}

compareToFileDistribution()

void compareToFileDistribution	(	char *	output,
		long	testCode,
		SDDS_DATASET *	SDDSin,
		char **	columnName,
		long	columnNames,
		char *	distFile,
		char *	distFileIndep,
		char *	distFileDepen )

Definition at line 314 of file sddsdistest.c.

                                                                                                                                                                                 {
  SDDS_DATASET SDDSdist;
  if (!SDDS_InitializeInput(&SDDSdist, distFile))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (SDDS_CheckColumn(&SDDSdist, distFileIndep, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY ||
      SDDS_CheckColumn(&SDDSdist, distFileDepen, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY)
    exit(EXIT_FAILURE);
  if (!SDDS_Terminate(&SDDSdist)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    exit(EXIT_FAILURE);
  }
  SDDS_Bomb("-fileDistribution option not implemented yet");
}

gaussianCDF()

double gaussianCDF

(

double

x

)

Definition at line 337 of file sddsdistest.c.

                             {
  long zSign;
  double z;
  z = (x - sampleMean) / sampleStDev;
  zSign = 1;
  if (z < 0) {
    zSign = -1;
    z = -z;
  }
  return (1 + zSign * erf(z / SQRT2)) / 2.0;
}

gaussianPDF()

double gaussianPDF

(

double

x

)

Definition at line 331 of file sddsdistest.c.

                             {
  double z;
  z = (x - sampleMean) / sampleStDev;
  return exp(-z * z / 2) / sqrt(PIx2);
}

ksTestWithFunction()

void ksTestWithFunction	(	double *	data,
		int64_t	rows,
		double(*	CDF )(double x),
		double *	statReturn,
		double *	sigLevelReturn )

Definition at line 536 of file sddsdistest.c.

                                                                                                                         {
  double CDF1, CDF2, dCDF1, dCDF2, CDF0, dCDFmaximum;
  int64_t row;
 
  qsort((void *)data, rows, sizeof(*data), double_cmpasc);
  dCDFmaximum = CDF1 = 0;
  for (row = 0; row < rows; row++) {
    CDF0 = (*CDF)(data[row]);
    CDF2 = (row + 1.0) / rows;
    dCDF1 = CDF0 - CDF1;
    if (dCDF1 < 0)
      dCDF1 = -dCDF1;
    dCDF2 = CDF0 - CDF2;
    if (dCDF2 < 0)
      dCDF2 = -dCDF2;
    CDF1 = CDF2;
    if (dCDF1 > dCDFmaximum)
      dCDFmaximum = dCDF1;
    if (dCDF2 > dCDFmaximum)
      dCDFmaximum = dCDF2;
  }
  *statReturn = dCDFmaximum;
  *sigLevelReturn = KS_Qfunction(sqrt((double)rows) * dCDFmaximum);
}

main()

int main	(	int	argc,
		char **	argv )

Definition at line 168 of file sddsdistest.c.

                                {
  int iArg;
  unsigned long dummyFlags, pipeFlags, majorOrderFlag;
  SCANNED_ARG *scanned;
  SDDS_DATASET SDDSin;
  char *input = NULL, *output = NULL, *distFile = NULL, **columnName = NULL, **sigmaName = NULL, **excludeName = NULL, *distFileIndep = NULL, *distFileDepen = NULL;
  long testCode = 0, distCode = -1, code, degreesFree = -1, columnNames = 0, excludeNames = 0;
  char *dofParameter = NULL;
  short columnMajorOrder = -1;
  int threads = 1;
 
  SDDS_RegisterProgramName(argv[0]);
  argc = scanargs(&scanned, argc, argv);
  if (argc < 3)
    bomb(NULL, USAGE);
 
  pipeFlags = 0;
 
  for (iArg = 1; iArg < argc; iArg++) {
    if (scanned[iArg].arg_type == OPTION) {
      /* process options here */
      code = match_string(scanned[iArg].list[0], option, N_OPTIONS, 0);
      switch (code) {
      case CLO_MAJOR_ORDER:
        majorOrderFlag = 0;
        scanned[iArg].n_items--;
        if (scanned[iArg].n_items > 0 &&
            !scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                          "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
                          "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)) {
          SDDS_Bomb("invalid -majorOrder syntax/values");
        }
        if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
          columnMajorOrder = 1;
        else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
          columnMajorOrder = 0;
        break;
      case CLO_PIPE:
        if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
          SDDS_Bomb("invalid -pipe syntax/values");
        break;
      case CLO_COLUMN:
        if ((scanned[iArg].n_items != 2 && scanned[iArg].n_items != 3) ||
            SDDS_StringIsBlank(scanned[iArg].list[1])) {
          SDDS_Bomb("invalid -column syntax/values");
        }
        moveToStringArray(&columnName, &columnNames, scanned[iArg].list + 1, 1);
        sigmaName = SDDS_Realloc(sigmaName, sizeof(*sigmaName) * columnNames);
        if (scanned[iArg].n_items == 3) {
          scanned[iArg].n_items -= 2;
          if (!scan_item_list(&dummyFlags, scanned[iArg].list + 2, &scanned[iArg].n_items, "sigma",
                              SDDS_STRING, sigmaName + columnNames - 1, 1, 1, NULL) ||
              dummyFlags != 1 ||
              SDDS_StringIsBlank(sigmaName[columnNames - 1])) {
            SDDS_Bomb("invalid -column syntax/values");
          }
        } else {
          sigmaName[columnNames - 1] = NULL;
        }
        break;
      case CLO_TEST:
        if (scanned[iArg].n_items != 2 ||
            (testCode = match_string(scanned[iArg].list[1], testChoice, N_TESTS, 0)) < 0) {
          SDDS_Bomb("invalid -test syntax/values");
        }
        break;
      case CLO_FILEDIST:
        if (scanned[iArg].n_items != 4) {
          SDDS_Bomb("too few qualifiers for -fileDistribution");
        }
        if (SDDS_StringIsBlank(distFile = scanned[iArg].list[1]) ||
            SDDS_StringIsBlank(distFileIndep = scanned[iArg].list[2]) ||
            SDDS_StringIsBlank(distFileDepen = scanned[iArg].list[3])) {
          SDDS_Bomb("invalid -fileDistribution values");
        }
        break;
      case CLO_GAUSSIAN:
      case CLO_POISSON:
      case CLO_STUDENT:
      case CLO_CHISQUARED:
        distCode = code;
        break;
      case CLO_DOF:
        if (scanned[iArg].n_items != 2) {
          SDDS_Bomb("too few qualifiers for -degreesOfFreedom");
        }
        if (scanned[iArg].list[1][0] == '@') {
          if (!SDDS_CopyString(&dofParameter, scanned[iArg].list[1] + 1)) {
            SDDS_Bomb("memory allocation failure");
          }
        } else if (sscanf(scanned[iArg].list[1], "%ld", &degreesFree) != 1 || degreesFree <= 1) {
          SDDS_Bomb("invalid degrees-of-freedom given for -student/-chiSquared");
        }
        break;
      case CLO_EXCLUDE:
        if (scanned[iArg].n_items < 2 || SDDS_StringIsBlank(scanned[iArg].list[1])) {
          SDDS_Bomb("invalid -exclude syntax/values");
        }
        moveToStringArray(&excludeName, &excludeNames, scanned[iArg].list + 1, scanned[iArg].n_items - 1);
        break;
      case CLO_THREADS:
        if (scanned[iArg].n_items != 2 ||
            !sscanf(scanned[iArg].list[1], "%d", &threads) ||
            threads < 1) {
          SDDS_Bomb("invalid -threads syntax");
        }
        break;
      default:
        fprintf(stderr, "error: unknown/ambiguous option: %s (%s)\n", scanned[iArg].list[0], argv[0]);
        exit(EXIT_FAILURE);
      }
    } else {
      if (!input)
        input = scanned[iArg].list[0];
      else if (!output)
        output = scanned[iArg].list[0];
      else
        SDDS_Bomb("too many filenames seen");
    }
  }
 
  processFilenames("sddsdistest", &input, &output, pipeFlags, 0, NULL);
 
  if (!SDDS_InitializeInput(&SDDSin, input))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (!columnNames)
    SDDS_Bomb("-column option must be supplied");
  if (!(columnNames = expandColumnPairNames(&SDDSin, &columnName, &sigmaName, columnNames, excludeName, excludeNames, FIND_NUMERIC_TYPE, 0))) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    SDDS_Bomb("named columns nonexistent or nonnumerical");
  }
  if (dofParameter && SDDS_CheckParameter(&SDDSin, dofParameter, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY)
    SDDS_Bomb("degrees-of-freedom parameter not found");
 
  if (distFile)
    compareToFileDistribution(output, testCode, &SDDSin, columnName, columnNames, distFile, distFileIndep, distFileDepen);
  else
    compareToDistribution(output, testCode, &SDDSin, columnName, sigmaName, columnNames, distCode, degreesFree, dofParameter, columnMajorOrder, threads);
 
  if (!SDDS_Terminate(&SDDSin)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}

poissonCDF()

double poissonCDF

(

double

xd

)

Definition at line 359 of file sddsdistest.c.

                             {
  double CDF, term, accuracy;
  long x, n;
 
  if (xd < 0)
    xd = 0;
  x = xd;
  xd = x;
  term = 1;
  CDF = 1;
  accuracy = POISSON_ACCURACY / exp(-sampleMean);
  for (n = 1; n <= x; n++) {
    term *= sampleMean / n;
    CDF += term;
    if (term < accuracy)
      break;
  }
  return CDF * exp(-sampleMean);
}

poissonPDF()

double poissonPDF

(

double

xd

)

Definition at line 351 of file sddsdistest.c.

                             {
  long x;
 
  if ((x = xd) < 0)
    return 0;
  return exp(-sampleMean + x * log(sampleMean) - lgamma(x + 1.0));
}

studentCDF()

double studentCDF

(

double

t

)

Definition at line 383 of file sddsdistest.c.

                            {
  if (t > 0)
    return 1 - betaInc(DOF / 2.0, 0.5, DOF / (DOF + t * t)) / 2;
  else
    return betaInc(DOF / 2.0, 0.5, DOF / (DOF + t * t)) / 2;
}

studentPDF()

double studentPDF

(

double

t

)

Definition at line 379 of file sddsdistest.c.

                            {
  return exp(-0.5 * (DOF + 1) * log(1 + t * t / DOF) + lgamma((DOF + 1.0) / 2.0) - lgamma(DOF / 2.0)) / sqrt(PI * DOF);
}