sharedjackknife.cpp

   1 /*
   2  *  sharedjackknife.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 3/30/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedjackknife.h"
  11
  12 /***************************************************************************************
  13
  14 ***************************************************************************************/
  15 double SharedJackknife::simpson(vector<int> abunds, double numInd, int numBins){
  16         double denom = numInd*(numInd-1);
  17         double sum = 0;
  18         for(int i = 0; i < numBins; i++)
  19                 sum += (double)abunds[i]*((double)abunds[i]-1)/denom;
  20
  21         return sum;
  22 }
  23
  24 /*****************************************************************************************/
  25
  26 double* SharedJackknife::jackknife(){
  27         int numBins = groups.at(0)->getNumBins()-1;
  28         vector<int> cArray(numBins);
  29         for(int i = 0; i < numBins; i++)
  30                 cArray[i] = 0;
  31
  32         double numInd = 0;
  33         for(int i = 0; i < numGroups; i++)
  34                 for(int j = 0; j < numBins; j++) {
  35                         int curAbund = groups.at(i)->get(j+1).abundance;
  36                         cArray[j] += curAbund;
  37                         numInd += (double)curAbund;
  38                 }
  39
  40         double baseD = 1/simpson(cArray, numInd, numBins);
  41
  42         vector<double> pseudoVals(numBins);
  43         double jackknifeEstimate = 0;
  44         for(int i = 0; i < numGroups; i++) {
  45                 for(int j = 0; j < numBins-1; j++) {
  46                         int abundDiff = -groups.at(i)->get(j+1).abundance;
  47                         if(i > 0)
  48                                 abundDiff += groups.at(i-1)->get(j+1).abundance;
  49
  50                         cArray[j] += abundDiff;
  51                         numInd += abundDiff;
  52                 }
  53
  54                 double curD = 1/simpson(cArray, numInd, numBins);
  55                 pseudoVals[i] = (double)numGroups*(baseD - curD) + curD;
  56                 jackknifeEstimate += pseudoVals[i];
  57         }
  58         jackknifeEstimate /= (double)numGroups;
  59
  60         double variance = 0;
  61         for(int i = 0; i < numGroups; i++)
  62                 variance += pow(pseudoVals[i]-jackknifeEstimate, 2);
  63
  64         variance /= (double)numGroups*((double)numGroups-1);
  65         double stErr = sqrt(variance);
  66         TDTable table;
  67         double confLimit = 0;
  68         if(numGroups <= 30)
  69                 confLimit = table.getConfLimit(numGroups-1, 1);
  70         else
  71                 confLimit = 1.645;
  72
  73         confLimit *= stErr;
  74
  75         double* rdata = new double[3];
  76         rdata[0] = baseD;
  77         rdata[1] = jackknifeEstimate - confLimit;
  78         rdata[2] = jackknifeEstimate + confLimit;
  79
  80         return rdata;
  81 }
  82
  83 /************************************************************************************************
  84 ************************************************************************************************/
  85
  86 EstOutput SharedJackknife::getValues(vector<SharedRAbundVector*> vectorShared){ //Fix this for collect, mistake was that it was made with summary in mind.
  87         try {
  88                 SharedRAbundVector* shared1 = vectorShared[0];
  89                 SharedRAbundVector* shared2 = vectorShared[1];
  90                 if(numGroups == -1) {
  91                         numGroups = m->getNumGroups();
  92                 }
  93
  94                 if(callCount == numGroups*(numGroups-1)/2) {
  95                         currentCallDone = true;
  96                         callCount = 0;
  97                 }
  98                 callCount++;
  99
 100                 if(currentCallDone) {
 101                         groups.clear();
 102                         currentCallDone = false;
 103                 }
 104
 105                 if(groups.size() != numGroups) {
 106                         if(groups.size() == 0)
 107                                 groups.push_back(shared1);
 108                         groups.push_back(shared2);
 109                 }
 110
 111                 if(groups.size() == numGroups && callCount < numGroups) {
 112                         data.resize(3,0);
 113
 114                         double* rdata = jackknife();
 115                         data[0] = rdata[0];
 116                         data[1] = rdata[1];
 117                         data[2] = rdata[2];
 118
 119                         if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
 120                         if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 121                         if (isnan(data[2]) || isinf(data[0])) { data[2] = 0; }
 122
 123                         return data;
 124                 }
 125
 126                 data.resize(3,0);
 127                 data[0] = 0;
 128                 data[1] = 0;
 129                 data[2] = 0;
 130
 131                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
 132                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 133                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
 134                 return data;
 135         }
 136
 137         catch(exception& e) {
 138                 m->errorOut(e, "SharedJackknife", "getValues");
 139                 exit(1);
 140         }
 141 }
 142
 143 /***********************************************************************/
 144