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