sharedkstest.cpp

   1 /*
   2  *  kstest.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 3/6/09.
   6  *  Copyright 2009 __MyCompanyName__. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedkstest.h"
  11
  12 /***********************************************************************/
  13
  14 EstOutput KSTest::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
  15         try {
  16                 data.resize(2,0);
  17
  18                 //Must return shared1 and shared2 to original order at conclusion of kstest
  19                 vector <individual> initData1 = shared1->getData();
  20                 vector <individual> initData2 = shared2->getData();
  21                 shared1->sortD();
  22                 shared2->sortD();
  23
  24                 int numNZ1 = shared1->numNZ();
  25                 int numNZ2 = shared2->numNZ();
  26                 double numInd1 = (double)shared1->getNumSeqs();
  27                 double numInd2 = (double)shared2->getNumSeqs();
  28
  29                 double maxDiff = -1;
  30                 double sum1 = 0;
  31                 double sum2 = 0;
  32                 for(int i = 1; i < shared1->getNumBins(); i++)
  33                 {
  34                         sum1 += shared1->get(i).abundance;
  35                         sum2 += shared2->get(i).abundance;
  36                         double diff = fabs((double)sum1/numInd1 - (double)sum2/numInd2);
  37                         if(diff > maxDiff)
  38                                 maxDiff = diff;
  39                 }
  40
  41                 double DStatistic = maxDiff*numNZ1*numNZ2;
  42                 double a = pow((double)(numNZ1 + numNZ2)/(numNZ1*numNZ2),.5);
  43                 //double pVal = exp(-2*pow(maxDiff/a,2));
  44                 double critVal = 1.36*a*numNZ1*numNZ2;
  45
  46                 /*cout << "Kolmogorov-Smirnov 2-sample test:\n";
  47                 if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25.
  48                         cout << "P-Value = " << pVal << "\nP-Value is the probability that the data sets are significantly different.\n";
  49                 else
  50                 {
  51                         //cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n";
  52                         cout << "D-Statistic = " << DStatistic << "\n";
  53                         cout << "95% Confidence Critical Value = " << critVal << "\n";
  54                         cout << "If D-Statistic is greater than the critical value then the data sets are significantly different at the 95% confidence level.\n\n";
  55                 }*/
  56
  57                 shared1->setData(initData1);
  58                 shared2->setData(initData2);
  59
  60                 data[0] = DStatistic;
  61                 data[1] = critVal;
  62                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
  63                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
  64
  65                 return data;
  66         }
  67         catch(exception& e) {
  68                 cout << "Standard Error: " << e.what() << " has occurred in the KSTest class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  69                 exit(1);
  70         }
  71         catch(...) {
  72                 cout << "An unknown error has occurred in the KSTest class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  73                 exit(1);
  74         }
  75 }
  76
  77 /***********************************************************************/