X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=calculator.cpp;h=3aa6a57a3d34f5e840b1a35d243b7c2903708358;hp=c45ddbffae174303b07e47b94d5183a37339022a;hb=d1c97b8c04bb75faca1e76ffad60b37a4d789d3d;hpb=c196b6b4768ccb84955d773ff0f22e4994d1ba7b diff --git a/calculator.cpp b/calculator.cpp index c45ddbf..3aa6a57 100644 --- a/calculator.cpp +++ b/calculator.cpp @@ -9,857 +9,59 @@ #include "calculator.h" -/***********************************************************************/ -void VecCalc::printElements(vector vec){ - try { - for(int i = 0; i < vec.size(); i++) - cout << vec.at(i) << " "; - cout << "\n\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} - -/***********************************************************************/ - -void VecCalc::printElements(vector vec){ - try { - for(int i = 0; i < vec.size(); i++) - cout << vec.at(i) << " "; - cout << "\n\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -int VecCalc::findString(vector vec, string str){ - try { - for(int i = 0; i < vec.size(); i++) - if(str.compare(vec.at(i)) == 0) - return i; - return -1; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double VecCalc::mean(vector vec){ - return sumElements(vec)/vec.size(); -} -/***********************************************************************/ -double VecCalc::stError(vector vec){ - try { - double sum = 0; - double m = mean(vec); - for(int i = 0; i < vec.size(); i++) - sum += pow(vec.at(i)-m,2); - return pow(sum/(vec.size()*(vec.size()-1)), .5); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} /***********************************************************************/ int VecCalc::sumElements(vector vec){ - try { return sumElements(vec,0); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } } /***********************************************************************/ int VecCalc::sumElements(vector vec, int index){ - try { + int sum = 0; for(int i = index; i < vec.size(); i++) sum += vec.at(i); return sum; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } + } /***********************************************************************/ double VecCalc::sumElements(vector vec){ - try { double sum = 0; for(int i = 0; i < vec.size(); i++) sum += vec.at(i); return sum; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } + } /***********************************************************************/ double VecCalc::sumElements(vector vec, int index){ - try { + double sum = 0; for(int i = index; i < vec.size(); i++) sum += vec.at(i); return sum; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double VecCalc::findMax(vector vec){ - try { - double max = -1000000.0; - for(int i = 0; i < vec.size(); i++) - if(vec.at(i) > max) - max = vec.at(i); - return max; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } + } /***********************************************************************/ int VecCalc::numNZ(vector vec){ - try { + int numNZ = 0; for(int i = 0; i < vec.size(); i++) if(vec.at(i) != 0) numNZ++; return numNZ; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } + } /***********************************************************************/ double VecCalc::numNZ(vector vec){ - try { + double numNZ = 0; for(int i = 0; i < vec.size(); i++) if(vec.at(i) != 0) numNZ++; return numNZ; } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double VecCalc::numPos(vector vec){ - try { - double numPos = 0; - for(int i = 0 ; i < vec.size(); i++) - if(vec.at(i) > 0) - numPos++; - return numPos; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double VecCalc::findMaxDiff(vector vec1, vector vec2){ - try { - double maxDiff = -10000000000.0; - for(int i = 0; i < vec1.size(); i++) - if(fabs(vec1.at(i)-vec2.at(i)) > maxDiff) - maxDiff = fabs(vec1.at(i)-vec2.at(i)); - return maxDiff; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double VecCalc::findDStat(vector vec1, vector vec2, double num){ - try { - double mDiff = findMaxDiff(vec1, vec2); - return (mDiff+.5)/num; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::findQuartiles(vector vec){ - try { - vector quartiles; - double max = vec.at(vec.size()-1); - double r1 = max/4; - double r2 = max*3/4; - bool r1found = false; - bool r2found = false; - for(int i = 0; i < vec.size(); i++) - { - if(vec.at(i) > r1 && !r1found) - { - quartiles.push_back(i); - r1found = true; - } - if(vec.at(i) > r2 && !r2found) - { - quartiles.push_back(i); - i = vec.size(); - } - } - return quartiles; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::add(vector vec, double x){ - try { - vector newVec; - for(int i = 0; i < vec.size(); i++) - newVec.push_back(vec.at(i)+x); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::multiply(vector vec, double x){ - try { - vector newVec; - for(int i = 0; i < vec.size(); i++) - newVec.push_back(vec.at(i)*x); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::power(vector vec, double p){ - try { - vector newVec; - for(int i = 0; i < vec.size(); i++) - newVec.push_back(pow(vec.at(i), p)); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::addVecs(vector vec1, vector vec2) //Vectors must be the same size. -{ try { - vector newVec; - for(int i = 0; i < vec1.size(); i++) - newVec.push_back(vec1.at(i)+vec2.at(i)); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::multVecs(vector vec1, vector vec2) //Vectors must be the same size. -{ try { - vector newVec; - for(int i = 0; i < vec1.size(); i++) - newVec.push_back(vec1.at(i)*vec2.at(i)); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::remDup(vector vec){ - try { - vector newVec; - newVec.push_back(vec.at(0)); - for(int i = 1; i < vec.size(); i++) - if(vec.at(i) != vec.at(i-1)) - newVec.push_back(vec.at(i)); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::genCVec(vector vec1){ - try { - vector vec2; - double sum = 0; - for(int i = 0; i < vec1.size(); i++) - { - sum += vec1.at(i); - vec2.push_back(sum); - } - return vec2; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::genRelVec(vector vec){ - try { - vector newVec; - double sum = sumElements(vec); - for(int i = 0; i < vec.size(); i++) - newVec.push_back(vec.at(i)/sum); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::genDiffVec(vector vec1, vector vec2){ - try { - vector newVec; - for(int i = 0; i < vec1.size(); i++) - newVec.push_back(vec1.at(i)-vec2.at(i)); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::genCSVec(vector vec){ - try { - vector newVec; - double curInd = vec.at(vec.size()-1); - double sumSpec = 0; - double cSumSpec = 0; - for(int i = vec.size()-1; i >= 0; i--) - { - if(vec.at(i) == curInd) - sumSpec++; - else - { - cSumSpec += sumSpec; - newVec.push_back(cSumSpec); - sumSpec = 1; - curInd = vec.at(i); - } - if(i == 0) - newVec.push_back(cSumSpec + sumSpec); - } - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::genTotVec(vector > vec){ - try { - vector newVec = vec.at(0); - for(int i = 1; i < vec.size(); i++) - newVec = addVecs(newVec, vec.at(i)); - return newVec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::quicksort(vector vec){ - try { - sort(vec.rbegin(), vec.rend()); - return vec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector > VecCalc::gen2DVec(vector vec, int div, int type){ - try { - vector > vec2D; - int gap = vec.size()/div; - for(int i = 0; i < div; i++) - { - vector inVec; - for(int j = 0; j < gap; j++) - if(type == 0) - inVec.push_back(vec.at(j + i*gap)); //Rows - else - inVec.push_back(vec.at(i + j*div)); //Columns - vec2D.push_back(inVec); - } - return vec2D; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -vector VecCalc::getSData(char name[]){ - try { - vector vec; - string data; - ifstream file(name); - if(file.is_open()) - { - while(file >> data) - vec.push_back(data); - file.close(); - } - return vec; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the VecCalc class function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ - -/***********************************************************************/ -void BrokenStick::doBStick(vector vec)//vec = The data vector. -{ try { - VecCalc vecCalc; - vec = vecCalc.quicksort(vec); - double numInd = vecCalc.sumElements(vec); - double numSpec = vecCalc.numNZ(vec); - vector cObsVec = vecCalc.genCVec(vec); - vector cExpVec; - vec = vecCalc.power(vec, -1); - double sumExp = 0; - for(int i = 0; i < numSpec; i++) - { - double n = 0; - for(int j = i; j < numSpec; j++) - n += 1/(double)(j+1); - - sumExp += numInd/numSpec*n; - cExpVec.push_back(sumExp); - } - - //Statistical analysis - double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec); - double DStatistic = maxDiff/numInd; - - cout << "D-Statistic = " << DStatistic << "\n"; - if(vec.size() > 30) - cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BrokenStick class Function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BrokenStick class function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} - -/***********************************************************************/ -double kEq(double k, double spec) -{ - return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec)); -} -/***********************************************************************/ -/*void GeometricSeries::doGeomTest(vector vec)//vec = the data vector -{ try { - VecCalc vecCalc; - vec = vecCalc.quicksort(vec);//Sorts vec - double numInd = vecCalc.sumElements(vec);//numInd = The total number of individuals in the data set. - double numSpec = vecCalc.numNZ(vec);//numSpec = the number of nonzero elements in the data set. - double min = -1.0*vecCalc.findMax(vecCalc.multiply(vec, -1)); - double k = .5; - double step = .49999; - while(fabs(min - numInd*kEq(k, numSpec)) > .0001)//This uses a binary search to find the value of k. - { - if(numInd*kEq(k, numSpec) > min) - k += step; - else - k -= step; - step /= 2; - } - cout << "k = " << k << "\n"; - double cK = 1/(1-pow(1-k, numSpec)); - - vector cObsVec = vecCalc.genCVec(vec); - vector cExpVec; - double sumExp = 0; - for(int i = 0; i < vec.size(); i++) - { - sumExp += numInd*cK*k*pow(1-k, i); - cExpVec.push_back(sumExp); - } - double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec); - double DStatistic = maxDiff/numInd; - - //Statistical Analysis - cout << "D-Statistic = " << DStatistic << "\n"; - if(vec.size() > 30) - cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the GeometricSeries class Function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the GeometricSeries class function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -}*/ - -/***********************************************************************/ -void Jackknifing::doJK(vector vec, double cols)//vec = the data vector if the data table was read left-to-right, top-to-bottom. cols = The number of columns in the data table. -{ try { - VecCalc vecCalc; - SSBPDiversityIndices ssbp; - double rows = vec.size()/cols; - vector > species = vecCalc.gen2DVec(vec, rows, 0);//converts vec into a 2 dimensional vector - vector sumRows; - vector pVals; - - for(int i = 0; i < rows; i++) - sumRows.push_back(vecCalc.sumElements(species.at(i))); - double st = 1/ssbp.getSimp(sumRows);//The observed estimate using the Simpson Index. Can be changed to use other indexes of diversity. - for(int i = 0; i < cols; i++) - { - vector newVec; - for(int j = 0; j < rows; j++) - newVec.push_back(vecCalc.sumElements(species.at(j))-species.at(j).at(i)); - pVals.push_back(cols*st-((cols-1)/ssbp.getSimp(newVec))); - } - - double mean = vecCalc.mean(pVals); - double stErr = vecCalc.stError(pVals);//stErr = standard error - TDTable table; - double confidence = .95; - double confLimit; - cout << "Confidence Level (.8, .9, .95, .98, .99, .998, .999): "; - cin >> confidence; - double confLevels[] = {.80,.90,.95,.98,.99,.998,.999}; - for(int i = 0; i < 7; i++) - if(confidence == confLevels[i]) - { - confLimit = table.getConfLimit(cols-2,i); - i = 7; - } - - //Statistical Analysis - cout << "Lower limit = " << mean - confLimit*stErr << "\n"; - cout << "Upper limit = " << mean + confLimit*stErr << "\n"; - cout << "Observed estimate = " << st << "\n\n"; - cout << "Jackknifed estimate = " << mean << "\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the Jackknifing class Function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the Jackknifing class function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -void KS2SampleTest::doKSTest(vector abun1, vector abun2)//abun1 = 1st vector of abundances, abun2 = 2nd vector of abundances -{ try { - VecCalc vecCalc; - double numNZ1 = vecCalc.numNZ(abun1);//Number of nonzero elements in abun1 - double numNZ2 = vecCalc.numNZ(abun2);//Number of nonzero elements in abun2 - abun1 = vecCalc.quicksort(abun1);//Sorts abun1 - abun2 = vecCalc.quicksort(abun2);//Sorts abun2 - abun1 = vecCalc.genRelVec(abun1);//Generates the relative vector for abun1 - abun2 = vecCalc.genRelVec(abun2);//Generates the relative vector for abun2 - - abun1 = vecCalc.genCVec(abun1);//Generates the cumulative vector for abun1 - abun2 = vecCalc.genCVec(abun2);//Generates the cumulative vector for abun2 - - double maxDiff = vecCalc.findMaxDiff(abun1, abun2); - double DStatistic = maxDiff*numNZ1*numNZ2; - - cout << "Null Hypothesis = There is no difference.\n"; - cout << "D-Statistic = " << DStatistic << "\n"; - - double a = pow((numNZ1 + numNZ2)/(numNZ1*numNZ2),.5); - double pVal = exp(-2*pow(maxDiff/a,2)); - - if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25. - cout << "P-Value = " << pVal << "\n\n"; - else - { - cout << "D-Statistic must be higher than the critical value to reject the null hypothesis.\n" ; - cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n"; - cout << "95% Confidence Critical Value = " << 1.36*a*numNZ1*numNZ2 << "\n"; - } - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the KS2SampleTest class Function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the KS2SampleTest class function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} - -/***********************************************************************/ - -void QStatistic::doQStat(vector vec)//vec = The data vector. -{ try { - VecCalc vecCalc; - vector cVec = vecCalc.genCSVec(vec); - vector q = vecCalc.findQuartiles(cVec);//q.at(0) is the index of the first quartile; q.at(1) is the index of the second quartile. - vector nDupVec = vecCalc.remDup(vec);//nDupVec only contains one of every unique element in cVec. - double Q; - if(q.at(0) != 0)//The case if neither quartile is 0 or 1 - Q = (.5*(cVec.at(q.at(0))-cVec.at(q.at(0)-1)) + (cVec.at(q.at(1)-1)-cVec.at(q.at(0))) + .5*(cVec.at(q.at(1))-cVec.at(q.at(1)-1)))/log(nDupVec.at(nDupVec.size()-1-q.at(1))/nDupVec.at(nDupVec.size()-1-q.at(0))); - else if(q.at(0) == 0 && (q.at(1) == 0 || q.at(1) == 1))//The case if the quartiles are both at index 0 or one is at 0 and the other at 1. - Q = (.5*cVec.at(0) + .5*(cVec.at(1)-cVec.at(0)))/log(nDupVec.at(nDupVec.size()-2)/nDupVec.at(nDupVec.size()-1)); - else if(q.at(0) == 0 && q.at(1) > 1) //The case if the lower quartile is at index 0 and upper quartile index is above index 1. - Q = (.5*cVec.at(0) + (cVec.at(q.at(1)-1)-cVec.at(q.at(0))) + .5*(cVec.at(q.at(1))-cVec.at(q.at(1)-1)))/log(nDupVec.at(nDupVec.size()-1-q.at(1))/nDupVec.at(nDupVec.size()-1)); - - cout << "Q = " << Q << "\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the QStatistic class Function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the QStatistic class function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double SSBPDiversityIndices::getShan(vector vec)//vec = The data vector. -{ try { - VecCalc vecCalc; - double nz = vecCalc.numNZ(vec); - double nSum = vecCalc.sumElements(vec); - double H = 0; - for(int i = 0; i < nz; i++) - H += vec.at(i)/nSum*log(vec.at(i)/nSum); - H *= -1; - return H; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the SSBPDiversityIndices class function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double SSBPDiversityIndices::getSimp(vector vec)//vec = The data vector. -{ try { - VecCalc vecCalc; - double nSum = vecCalc.sumElements(vec); - double D = 0; - for(int j = 0; j < vec.size(); j++) - D += vec.at(j)*(vec.at(j)-1)/(nSum*(nSum-1)); - return D; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the SSBPDiversityIndices class function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double SSBPDiversityIndices::getBP(vector vec)//vec = The data vector. -{ try { - VecCalc vecCalc; - double nSum = vecCalc.sumElements(vec); - return vecCalc.findMax(vec)/nSum; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the SSBPDiversityIndices class function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -void SSBPDiversityIndices::doSSBP(vector vec)//vec = The data vector. -{ try { - VecCalc vecCalc; - double nz = vecCalc.numNZ(vec); - - //Shannon index - double H = getShan(vec); - cout << "H = " << H << "\n"; - cout << "Eveness = " << H/log(nz) << "\n\n"; - - //Simpson index - double D = getSimp(vec); - cout << "D diversity = " << 1/D << "\n"; - cout << "Eveness = " << 1/D/nz << "\n\n"; - - //Berger-Parker index - double BP = getBP(vec); - cout << "BP index = " << BP << "\n\n"; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the SSBPDiversityIndices class function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} /***********************************************************************/ double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom -{ try { +{ //Found on http://www.vgtu.lt/leidiniai/elektroniniai/Probability.pdf/Table%203.pdf //Confidence Level .90 .95 .975 .99 .995 .999 .9995 @@ -895,124 +97,10 @@ double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646}}; return values[row][col]; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } + } /***********************************************************************/ -double KOSTable::getConfLimit(int index) //Table of Critical values for N = 4-20 for One-Sample Kolmogorov-Smirnov Test -{ try { - //Confidence Level = .05 - //Sample size = 4-20. - //Found on http://www.utdallas.edu/~herve/MolinAbdi1998-LillieforsTechReport.pdf - double values[21] = {.9011, - .6372, - .5203, - .4506, - 0.3754, - 0.3427, - 0.3245, - 0.3041, - 0.2875, - 0.2744, - 0.2616, - 0.2506, - 0.2426, - 0.2337, - 0.2257, - 0.2196, - 0.2128, - 0.2071, - 0.2018, - 0.1965, - 0.1920, - }; - return values[index]; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/*********************************************************************** -void TrunLN::doTrunLN(vector indVec, vector specVec) //indVec = individuals vector, specVec = species vector -{ - VecCalc vecCalc; - double numSpec = vecCalc.sumElements(specVec); //numSpec = The total number of species - double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec,specVec)); //numInd = The total number of individuals - - double obsMean = 0; - for(int i = 0; i < indVec.size(); i++) - obsMean += log10(indVec.at(i)); - obsMean /= numSpec; //obsMean = observed mean of the individuals vector - cout << "obsMean = " << obsMean << "\n"; - double variance = 0; - for(int t = 0; t < indVec.size(); t++) - variance += pow(log10(indVec.at(t))-obsMean,2)/numSpec; - - double rO = 0; - for(int k = 0; k < indVec.size(); k++) - rO += log10(indVec.at(k)); - rO /= indVec.size(); - double veilLine = .5;//The desired veil line. - double auxFunc = -(obsMean-rO)/(obsMean-log10(veilLine)); - double uX = obsMean-auxFunc*(obsMean-log10(veilLine)); - double vX = variance + auxFunc*pow(obsMean-log10(veilLine),2); - double z = (log10(veilLine)-uX)/pow(vX, .5); - double p = .5*(erf(z)+1); - double specRichness = numSpec/(1-p); - - double numUnseen = .5*(erf((log10(.5)-uX)/pow(vX,.5))+1)*specRichness; - - - vector cExp; - for(int i = 1; i < 8; i++) - { - double a = pow(10, i-1)+.5; - double b = log10(a); - double c = (b - uX)/pow(vX,.5); - double d = .5*(erf(c)+1); - double numS = d*specRichness; - double toPush = numS - numUnseen; - cExp.push_back(toPush); - } - vector cObs; - double sumOct = 0; - for(int i = 0; i < 8; i++) - { - sumOct = 0; - for(int r = 0; r < indVec.size(); r++) - { - if(indVec.at(r) < pow(10, i-1)+.5) - sumOct += specVec.at(r); - else - { - cObs.push_back(sumOct); - sumOct = specVec.at(r); - r = indVec.size(); - } - if(r == indVec.size()-1) - cObs.push_back(sumOct); - } - } - //Statistical Analysis - double d = vecCalc.findDStat(cExp, cObs, numSpec); - cout << "DStat = " << d << "\n"; - cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n"; - cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n"; -} -/***********************************************************************/