geom.cpp

   1 /*
   2  *  geom.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 2/23/09.
   6  *  Copyright 2009 __MyCompanyName__. All rights reserved.
   7  *
   8  */
   9
  10 #include "geom.h"
  11 #include "calculator.h"
  12
  13 /***********************************************************************/
  14
  15 double Geom::kEq(double k, double spec){
  16         return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec));
  17 }
  18
  19 RAbundVector Geom::getRAbundVector(SAbundVector* rank){
  20                 vector <int> rData;
  21                 int mr = 1;
  22                 int nb = 0;
  23                 int ns = 0;
  24
  25                 for(int i = rank->size()-1; i > 0; i--)
  26                 {
  27                         int cur = rank->get(i);
  28                         if(mr == 1 && cur > 0)
  29                                 mr = i;
  30                         nb += cur;
  31                         ns += i*cur;
  32                         for(int j = 0; j < cur; j++)
  33                                 rData.push_back(i);
  34                 }
  35
  36                 RAbundVector rav = RAbundVector(rData, mr, nb, ns);
  37                 rav.setLabel(rank->getLabel());
  38                 return rav;
  39 }
  40
  41 /***********************************************************************************/
  42
  43 /***********************************************************************************/
  44 EstOutput Geom::getValues(SAbundVector* rank){
  45         try {
  46                 data.resize(2,0);
  47
  48                 rdata = getRAbundVector(rank);
  49                 int numInd = rdata.getNumSeqs();
  50                 int numSpec = rdata.getNumBins();
  51                 int min = rdata.get(rdata.size()-1);
  52                 double k = .5;
  53                 double step = .49999;
  54
  55                 while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) //This uses a binary search to find the value of k.
  56                 {
  57                         if(numInd*kEq(k, numSpec) > min)
  58                                 k += step;
  59                         else
  60                                 k -= step;
  61                         step /= 2;
  62                 }
  63
  64                 double cK = 1/(1-pow(1-k, numSpec));
  65                 double sumExp = 0;
  66                 double sumObs = 0;
  67                 double maxDiff = 0;
  68
  69                 for(int i = 0; i < rdata.size(); i++)
  70                 {
  71                         sumObs += rdata.get(i);
  72                         sumExp += numInd*cK*k*pow(1-k, i);
  73                         double diff = fabs(sumObs-sumExp);
  74                         if(diff > maxDiff)
  75                                 maxDiff = diff;
  76                 }
  77
  78                 double DStatistic = maxDiff/numInd;
  79                 double critVal = 0;
  80                 /*cout << "Geom:\n";
  81                 cout << "D-Statistic = " << DStatistic << "\n";
  82                 cout << "Critical value for 95% confidence interval = ";*/
  83                 if(rdata.size() > 20)
  84                 {
  85                         critVal = .886/sqrt(rdata.size());
  86                 }
  87                 else
  88                 {
  89                         KOSTable table;
  90                         critVal = table.getConfLimit(numSpec);
  91                 }
  92                 /*cout << critVal << "\n";
  93                 cout << "If D-Statistic is less than the critical value then the data fits the Geometric Series model w/ 95% confidence.\n\n";*/
  94
  95                 data[0] = DStatistic;
  96                 data[1] = critVal;
  97
  98                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
  99                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 100
 101                 return data;
 102         }
 103         catch(exception& e) {
 104                 cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 105                 exit(1);
 106         }
 107         catch(...) {
 108                 cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 109                 exit(1);
 110         }
 111 }
 112
 113 /***********************************************************************/
 114
 115
 116
 117
 118
 119
 120