sharedkstest.cpp

   1 /*
   2  *  kstest.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 3/6/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedkstest.h"
  11
  12 /***********************************************************************/
  13
  14 EstOutput KSTest::getValues(vector<SharedRAbundVector*> shared){
  15         try {
  16                 data.resize(3,0);
  17
  18                 //Must return shared1 and shared2 to original order at conclusion of kstest
  19                 vector <individual> initData1 = shared[0]->getData();
  20                 vector <individual> initData2 = shared[1]->getData();
  21                 shared[0]->sortD();
  22                 shared[1]->sortD();
  23
  24                 int numNZ1 = shared[0]->numNZ();
  25                 int numNZ2 = shared[1]->numNZ();
  26                 double numInd1 = (double)shared[0]->getNumSeqs();
  27                 double numInd2 = (double)shared[1]->getNumSeqs();
  28
  29                 double maxDiff = -1;
  30                 double sum1 = 0;
  31                 double sum2 = 0;
  32                 for(int i = 1; i < shared[0]->getNumBins(); i++)
  33                 {
  34                         sum1 += shared[0]->get(i).abundance;
  35                         sum2 += shared[1]->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                 shared[0]->setData(initData1);
  47                 shared[1]->setData(initData2);
  48
  49                 data[0] = DStatistic;
  50                 data[1] = critVal;
  51                 data[2] = 0;
  52
  53                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
  54                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
  55
  56                 return data;
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "KSTest", "getValues");
  60                 exit(1);
  61         }
  62 }
  63
  64 /***********************************************************************/