5 * Created by Sarah Westcott on 11/18/08.
6 * Copyright 2008 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "calculator.h"
12 /***********************************************************************
13 void VecCalc::printElements(vector<double> vec){
15 for(int i = 0; i < vec.size(); i++)
16 cout << vec.at(i) << " ";
20 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
24 cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
29 /***********************************************************************
31 void VecCalc::printElements(vector<string> vec){
33 for(int i = 0; i < vec.size(); i++)
34 cout << vec.at(i) << " ";
38 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
42 cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
46 /***********************************************************************
47 int VecCalc::findString(vector<string> vec, string str){
49 for(int i = 0; i < vec.size(); i++)
50 if(str.compare(vec.at(i)) == 0)
55 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
59 cout << "An unknown error has occurred in the VecCalc class function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
63 /***********************************************************************
64 double VecCalc::mean(vector<double> vec){
65 return sumElements(vec)/vec.size();
67 /***********************************************************************
68 double VecCalc::stError(vector<double> vec){
72 for(int i = 0; i < vec.size(); i++)
73 sum += pow(vec.at(i)-m,2);
74 return pow(sum/(vec.size()*(vec.size()-1)), .5);
77 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
81 cout << "An unknown error has occurred in the VecCalc class function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
85 /***********************************************************************/
86 int VecCalc::sumElements(vector<int> vec){
88 return sumElements(vec,0);
91 errorOut(e, "VecCalc", "sumElements");
95 /***********************************************************************/
96 int VecCalc::sumElements(vector<int> vec, int index){
99 for(int i = index; i < vec.size(); i++)
103 catch(exception& e) {
104 errorOut(e, "VecCalc", "sumElements");
109 /***********************************************************************/
110 double VecCalc::sumElements(vector<double> vec){
113 for(int i = 0; i < vec.size(); i++)
117 catch(exception& e) {
118 errorOut(e, "VecCalc", "sumElements");
122 /***********************************************************************/
123 double VecCalc::sumElements(vector<double> vec, int index){
126 for(int i = index; i < vec.size(); i++)
130 catch(exception& e) {
131 errorOut(e, "VecCalc", "sumElements");
135 /***********************************************************************
136 double VecCalc::findMax(vector<double> vec){
138 double max = -1000000.0;
139 for(int i = 0; i < vec.size(); i++)
144 catch(exception& e) {
145 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
149 cout << "An unknown error has occurred in the VecCalc class function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
153 /***********************************************************************/
154 int VecCalc::numNZ(vector<int> vec){
157 for(int i = 0; i < vec.size(); i++)
162 catch(exception& e) {
163 errorOut(e, "VecCalc", "numNZ");
167 /***********************************************************************/
168 double VecCalc::numNZ(vector<double> vec){
171 for(int i = 0; i < vec.size(); i++)
176 catch(exception& e) {
177 errorOut(e, "VecCalc", "numNZ");
181 /***********************************************************************
182 double VecCalc::numPos(vector<double> vec){
185 for(int i = 0 ; i < vec.size(); i++)
190 catch(exception& e) {
191 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
195 cout << "An unknown error has occurred in the VecCalc class function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
199 /***********************************************************************
200 double VecCalc::findMaxDiff(vector<double> vec1, vector<double> vec2){
202 double maxDiff = -10000000000.0;
203 for(int i = 0; i < vec1.size(); i++)
204 if(fabs(vec1.at(i)-vec2.at(i)) > maxDiff)
205 maxDiff = fabs(vec1.at(i)-vec2.at(i));
208 catch(exception& e) {
209 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
213 cout << "An unknown error has occurred in the VecCalc class function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
217 /***********************************************************************
218 double VecCalc::findDStat(vector<double> vec1, vector<double> vec2, double num){
220 double mDiff = findMaxDiff(vec1, vec2);
221 return (mDiff+.5)/num;
223 catch(exception& e) {
224 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
228 cout << "An unknown error has occurred in the VecCalc class function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
232 /***********************************************************************
233 vector<int> VecCalc::findQuartiles(vector<double> vec){
235 vector<int> quartiles;
236 double max = vec.at(vec.size()-1);
239 bool r1found = false;
240 bool r2found = false;
241 for(int i = 0; i < vec.size(); i++)
243 if(vec.at(i) > r1 && !r1found)
245 quartiles.push_back(i);
248 if(vec.at(i) > r2 && !r2found)
250 quartiles.push_back(i);
256 catch(exception& e) {
257 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
261 cout << "An unknown error has occurred in the VecCalc class function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
265 /***********************************************************************
266 vector<double> VecCalc::add(vector<double> vec, double x){
268 vector<double> newVec;
269 for(int i = 0; i < vec.size(); i++)
270 newVec.push_back(vec.at(i)+x);
273 catch(exception& e) {
274 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
278 cout << "An unknown error has occurred in the VecCalc class function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
282 /***********************************************************************
283 vector<double> VecCalc::multiply(vector<double> vec, double x){
285 vector<double> newVec;
286 for(int i = 0; i < vec.size(); i++)
287 newVec.push_back(vec.at(i)*x);
290 catch(exception& e) {
291 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
295 cout << "An unknown error has occurred in the VecCalc class function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
299 /***********************************************************************
300 vector<double> VecCalc::power(vector<double> vec, double p){
302 vector<double> newVec;
303 for(int i = 0; i < vec.size(); i++)
304 newVec.push_back(pow(vec.at(i), p));
307 catch(exception& e) {
308 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
312 cout << "An unknown error has occurred in the VecCalc class function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
316 /***********************************************************************
317 vector<double> VecCalc::addVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
319 vector<double> newVec;
320 for(int i = 0; i < vec1.size(); i++)
321 newVec.push_back(vec1.at(i)+vec2.at(i));
324 catch(exception& e) {
325 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
329 cout << "An unknown error has occurred in the VecCalc class function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
333 /***********************************************************************
334 vector<double> VecCalc::multVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
336 vector<double> newVec;
337 for(int i = 0; i < vec1.size(); i++)
338 newVec.push_back(vec1.at(i)*vec2.at(i));
341 catch(exception& e) {
342 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
346 cout << "An unknown error has occurred in the VecCalc class function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
350 /***********************************************************************
351 vector<double> VecCalc::remDup(vector<double> vec){
353 vector<double> newVec;
354 newVec.push_back(vec.at(0));
355 for(int i = 1; i < vec.size(); i++)
356 if(vec.at(i) != vec.at(i-1))
357 newVec.push_back(vec.at(i));
360 catch(exception& e) {
361 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
365 cout << "An unknown error has occurred in the VecCalc class function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
369 /***********************************************************************
370 vector<double> VecCalc::genCVec(vector<double> vec1){
374 for(int i = 0; i < vec1.size(); i++)
381 catch(exception& e) {
382 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
386 cout << "An unknown error has occurred in the VecCalc class function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
390 /***********************************************************************
391 vector<double> VecCalc::genRelVec(vector<double> vec){
393 vector<double> newVec;
394 double sum = sumElements(vec);
395 for(int i = 0; i < vec.size(); i++)
396 newVec.push_back(vec.at(i)/sum);
399 catch(exception& e) {
400 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
404 cout << "An unknown error has occurred in the VecCalc class function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
408 /***********************************************************************
409 vector<double> VecCalc::genDiffVec(vector<double> vec1, vector<double> vec2){
411 vector<double> newVec;
412 for(int i = 0; i < vec1.size(); i++)
413 newVec.push_back(vec1.at(i)-vec2.at(i));
416 catch(exception& e) {
417 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
421 cout << "An unknown error has occurred in the VecCalc class function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
425 /***********************************************************************
426 vector<double> VecCalc::genCSVec(vector<double> vec){
428 vector<double> newVec;
429 double curInd = vec.at(vec.size()-1);
432 for(int i = vec.size()-1; i >= 0; i--)
434 if(vec.at(i) == curInd)
439 newVec.push_back(cSumSpec);
444 newVec.push_back(cSumSpec + sumSpec);
448 catch(exception& e) {
449 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
453 cout << "An unknown error has occurred in the VecCalc class function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
457 /***********************************************************************
458 vector<double> VecCalc::genTotVec(vector<vector<double> > vec){
460 vector<double> newVec = vec.at(0);
461 for(int i = 1; i < vec.size(); i++)
462 newVec = addVecs(newVec, vec.at(i));
465 catch(exception& e) {
466 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
470 cout << "An unknown error has occurred in the VecCalc class function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
474 /***********************************************************************
475 vector<double> VecCalc::quicksort(vector<double> vec){
477 sort(vec.rbegin(), vec.rend());
480 catch(exception& e) {
481 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
485 cout << "An unknown error has occurred in the VecCalc class function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
489 /***********************************************************************
490 vector<vector<double> > VecCalc::gen2DVec(vector<double> vec, int div, int type){
492 vector<vector<double> > vec2D;
493 int gap = vec.size()/div;
494 for(int i = 0; i < div; i++)
496 vector<double> inVec;
497 for(int j = 0; j < gap; j++)
499 inVec.push_back(vec.at(j + i*gap)); //Rows
501 inVec.push_back(vec.at(i + j*div)); //Columns
502 vec2D.push_back(inVec);
506 catch(exception& e) {
507 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
511 cout << "An unknown error has occurred in the VecCalc class function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
515 /***********************************************************************
516 vector<string> VecCalc::getSData(char name[]){
529 catch(exception& e) {
530 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
534 cout << "An unknown error has occurred in the VecCalc class function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
538 /***********************************************************************/
540 /***********************************************************************
541 void BrokenStick::doBStick(vector<double> vec)//vec = The data vector.
544 vec = vecCalc.quicksort(vec);
545 double numInd = vecCalc.sumElements(vec);
546 double numSpec = vecCalc.numNZ(vec);
547 vector<double> cObsVec = vecCalc.genCVec(vec);
548 vector<double> cExpVec;
549 vec = vecCalc.power(vec, -1);
551 for(int i = 0; i < numSpec; i++)
554 for(int j = i; j < numSpec; j++)
555 n += 1/(double)(j+1);
557 sumExp += numInd/numSpec*n;
558 cExpVec.push_back(sumExp);
561 //Statistical analysis
562 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
563 double DStatistic = maxDiff/numInd;
565 cout << "D-Statistic = " << DStatistic << "\n";
567 cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
569 catch(exception& e) {
570 cout << "Standard Error: " << e.what() << " has occurred in the BrokenStick class Function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
574 cout << "An unknown error has occurred in the BrokenStick class function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
579 /***********************************************************************
580 double kEq(double k, double spec)
582 return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec));
584 /***********************************************************************/
585 /*void GeometricSeries::doGeomTest(vector<double> vec)//vec = the data vector
588 vec = vecCalc.quicksort(vec);//Sorts vec
589 double numInd = vecCalc.sumElements(vec);//numInd = The total number of individuals in the data set.
590 double numSpec = vecCalc.numNZ(vec);//numSpec = the number of nonzero elements in the data set.
591 double min = -1.0*vecCalc.findMax(vecCalc.multiply(vec, -1));
593 double step = .49999;
594 while(fabs(min - numInd*kEq(k, numSpec)) > .0001)//This uses a binary search to find the value of k.
596 if(numInd*kEq(k, numSpec) > min)
602 cout << "k = " << k << "\n";
603 double cK = 1/(1-pow(1-k, numSpec));
605 vector<double> cObsVec = vecCalc.genCVec(vec);
606 vector<double> cExpVec;
608 for(int i = 0; i < vec.size(); i++)
610 sumExp += numInd*cK*k*pow(1-k, i);
611 cExpVec.push_back(sumExp);
613 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
614 double DStatistic = maxDiff/numInd;
616 //Statistical Analysis
617 cout << "D-Statistic = " << DStatistic << "\n";
619 cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
621 catch(exception& e) {
622 cout << "Standard Error: " << e.what() << " has occurred in the GeometricSeries class Function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
626 cout << "An unknown error has occurred in the GeometricSeries class function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
631 /***********************************************************************
632 void Jackknifing::doJK(vector<double> 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.
635 SSBPDiversityIndices ssbp;
636 double rows = vec.size()/cols;
637 vector<vector<double> > species = vecCalc.gen2DVec(vec, rows, 0);//converts vec into a 2 dimensional vector
638 vector<double> sumRows;
639 vector<double> pVals;
641 for(int i = 0; i < rows; i++)
642 sumRows.push_back(vecCalc.sumElements(species.at(i)));
643 double st = 1/ssbp.getSimp(sumRows);//The observed estimate using the Simpson Index. Can be changed to use other indexes of diversity.
644 for(int i = 0; i < cols; i++)
646 vector<double> newVec;
647 for(int j = 0; j < rows; j++)
648 newVec.push_back(vecCalc.sumElements(species.at(j))-species.at(j).at(i));
649 pVals.push_back(cols*st-((cols-1)/ssbp.getSimp(newVec)));
652 double mean = vecCalc.mean(pVals);
653 double stErr = vecCalc.stError(pVals);//stErr = standard error
655 double confidence = .95;
657 cout << "Confidence Level (.8, .9, .95, .98, .99, .998, .999): ";
659 double confLevels[] = {.80,.90,.95,.98,.99,.998,.999};
660 for(int i = 0; i < 7; i++)
661 if(confidence == confLevels[i])
663 confLimit = table.getConfLimit(cols-2,i);
667 //Statistical Analysis
668 cout << "Lower limit = " << mean - confLimit*stErr << "\n";
669 cout << "Upper limit = " << mean + confLimit*stErr << "\n";
670 cout << "Observed estimate = " << st << "\n\n";
671 cout << "Jackknifed estimate = " << mean << "\n";
673 catch(exception& e) {
674 cout << "Standard Error: " << e.what() << " has occurred in the Jackknifing class Function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
678 cout << "An unknown error has occurred in the Jackknifing class function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
682 /***********************************************************************
683 void KS2SampleTest::doKSTest(vector<double> abun1, vector<double> abun2)//abun1 = 1st vector of abundances, abun2 = 2nd vector of abundances
686 double numNZ1 = vecCalc.numNZ(abun1);//Number of nonzero elements in abun1
687 double numNZ2 = vecCalc.numNZ(abun2);//Number of nonzero elements in abun2
688 abun1 = vecCalc.quicksort(abun1);//Sorts abun1
689 abun2 = vecCalc.quicksort(abun2);//Sorts abun2
690 abun1 = vecCalc.genRelVec(abun1);//Generates the relative vector for abun1
691 abun2 = vecCalc.genRelVec(abun2);//Generates the relative vector for abun2
693 abun1 = vecCalc.genCVec(abun1);//Generates the cumulative vector for abun1
694 abun2 = vecCalc.genCVec(abun2);//Generates the cumulative vector for abun2
696 double maxDiff = vecCalc.findMaxDiff(abun1, abun2);
697 double DStatistic = maxDiff*numNZ1*numNZ2;
699 cout << "Null Hypothesis = There is no difference.\n";
700 cout << "D-Statistic = " << DStatistic << "\n";
702 double a = pow((numNZ1 + numNZ2)/(numNZ1*numNZ2),.5);
703 double pVal = exp(-2*pow(maxDiff/a,2));
705 if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25.
706 cout << "P-Value = " << pVal << "\n\n";
709 cout << "D-Statistic must be higher than the critical value to reject the null hypothesis.\n" ;
710 cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n";
711 cout << "95% Confidence Critical Value = " << 1.36*a*numNZ1*numNZ2 << "\n";
714 catch(exception& e) {
715 cout << "Standard Error: " << e.what() << " has occurred in the KS2SampleTest class Function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
719 cout << "An unknown error has occurred in the KS2SampleTest class function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
724 /***********************************************************************
726 void QStatistic::doQStat(vector<double> vec)//vec = The data vector.
729 vector<double> cVec = vecCalc.genCSVec(vec);
730 vector<int> q = vecCalc.findQuartiles(cVec);//q.at(0) is the index of the first quartile; q.at(1) is the index of the second quartile.
731 vector<double> nDupVec = vecCalc.remDup(vec);//nDupVec only contains one of every unique element in cVec.
733 if(q.at(0) != 0)//The case if neither quartile is 0 or 1
734 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)));
735 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.
736 Q = (.5*cVec.at(0) + .5*(cVec.at(1)-cVec.at(0)))/log(nDupVec.at(nDupVec.size()-2)/nDupVec.at(nDupVec.size()-1));
737 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.
738 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));
740 cout << "Q = " << Q << "\n";
742 catch(exception& e) {
743 cout << "Standard Error: " << e.what() << " has occurred in the QStatistic class Function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
747 cout << "An unknown error has occurred in the QStatistic class function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
751 /***********************************************************************
752 double SSBPDiversityIndices::getShan(vector<double> vec)//vec = The data vector.
755 double nz = vecCalc.numNZ(vec);
756 double nSum = vecCalc.sumElements(vec);
758 for(int i = 0; i < nz; i++)
759 H += vec.at(i)/nSum*log(vec.at(i)/nSum);
763 catch(exception& e) {
764 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
768 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
772 /***********************************************************************
773 double SSBPDiversityIndices::getSimp(vector<double> vec)//vec = The data vector.
776 double nSum = vecCalc.sumElements(vec);
778 for(int j = 0; j < vec.size(); j++)
779 D += vec.at(j)*(vec.at(j)-1)/(nSum*(nSum-1));
782 catch(exception& e) {
783 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
787 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
791 /***********************************************************************
792 double SSBPDiversityIndices::getBP(vector<double> vec)//vec = The data vector.
795 double nSum = vecCalc.sumElements(vec);
796 return vecCalc.findMax(vec)/nSum;
798 catch(exception& e) {
799 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
803 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
807 /***********************************************************************
808 void SSBPDiversityIndices::doSSBP(vector<double> vec)//vec = The data vector.
811 double nz = vecCalc.numNZ(vec);
814 double H = getShan(vec);
815 cout << "H = " << H << "\n";
816 cout << "Eveness = " << H/log(nz) << "\n\n";
819 double D = getSimp(vec);
820 cout << "D diversity = " << 1/D << "\n";
821 cout << "Eveness = " << 1/D/nz << "\n\n";
823 //Berger-Parker index
824 double BP = getBP(vec);
825 cout << "BP index = " << BP << "\n\n";
827 catch(exception& e) {
828 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
832 cout << "An unknown error has occurred in the SSBPDiversityIndices class function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
836 /***********************************************************************/
837 double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom
839 //Found on http://www.vgtu.lt/leidiniai/elektroniniai/Probability.pdf/Table%203.pdf
841 //Confidence Level .90 .95 .975 .99 .995 .999 .9995
842 double values[30][7] = {{3.078, 6.314, 12.706, 31.821, 63.656, 318.289, 636.578},
843 {1.886, 2.920, 4.303, 6.965, 9.925, 22.328, 31.600},
844 {1.638, 2.353, 3.182, 4.541, 5.841, 10.214, 12.924},
845 {1.533, 2.132, 2.776, 3.747, 4.604, 7.173, 8.610},
846 {1.476, 2.015, 2.571, 3.365, 4.032, 5.894, 6.869},
847 {1.440, 1.943, 2.447, 3.143, 3.707, 5.208, 5.959},
848 {1.415, 1.895, 2.365, 2.998, 3.499, 4.785, 5.408},
849 {1.397, 1.860, 2.306, 2.896, 3.355, 4.501, 5.041},
850 {1.383, 1.833, 2.262, 2.821, 3.250, 4.297, 4.781},
851 {1.372, 1.812, 2.228, 2.764, 3.169, 4.144, 4.587},
852 {1.363, 1.796, 2.201, 2.718, 3.106, 4.025, 4.437},
853 {1.356, 1.782, 2.179, 2.681, 3.055, 3.930, 4.318},
854 {1.350, 1.771, 2.160, 2.650, 3.012, 3.852, 4.221},
855 {1.345, 1.761, 2.145, 2.624, 2.977, 3.787, 4.140},
856 {1.341, 1.753, 2.131, 2.602, 2.947, 3.733, 4.073},
857 {1.337, 1.746, 2.120, 2.583, 2.921, 3.686, 4.015},
858 {1.333, 1.740, 2.110, 2.567, 2.898, 3.646, 3.965},
859 {1.330, 1.734, 2.101, 2.552, 2.878, 3.610, 3.922},
860 {1.328, 1.729, 2.093, 2.539, 2.861, 3.579, 3.883},
861 {1.325, 1.725, 2.086, 2.528, 2.845, 3.552, 3.850},
862 {1.323, 1.721, 2.080, 2.518, 2.831, 3.527, 3.819},
863 {1.321, 1.717, 2.074, 2.508, 2.819, 3.505, 3.792},
864 {1.319, 1.714, 2.069, 2.500, 2.807, 3.485, 3.768},
865 {1.318, 1.711, 2.064, 2.492, 2.797, 3.467, 3.745},
866 {1.316, 1.708, 2.060, 2.485, 2.787, 3.450, 3.725},
867 {1.315, 1.706, 2.056, 2.479, 2.779, 3.435, 3.707},
868 {1.314, 1.703, 2.052, 2.473, 2.771, 3.421, 3.689},
869 {1.313, 1.701, 2.048, 2.467, 2.763, 3.408, 3.674},
870 {1.311, 1.699, 2.045, 2.462, 2.756, 3.396, 3.660},
871 {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646}};
873 return values[row][col];
875 catch(exception& e) {
876 errorOut(e, "TDTable", "getConfLimit");
881 /***********************************************************************
882 double KOSTable::getConfLimit(int index) //Table of Critical values for N = 4-20 for One-Sample Kolmogorov-Smirnov Test
884 //Confidence Level = .05
885 //Sample size = 4-20.
886 //Found on http://www.utdallas.edu/~herve/MolinAbdi1998-LillieforsTechReport.pdf
888 double values[21] = {.9011,
910 return values[index];
912 catch(exception& e) {
913 cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
917 cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
922 /***********************************************************************
923 void TrunLN::doTrunLN(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
926 double numSpec = vecCalc.sumElements(specVec); //numSpec = The total number of species
927 double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec,specVec)); //numInd = The total number of individuals
930 for(int i = 0; i < indVec.size(); i++)
931 obsMean += log10(indVec.at(i));
932 obsMean /= numSpec; //obsMean = observed mean of the individuals vector
933 cout << "obsMean = " << obsMean << "\n";
935 for(int t = 0; t < indVec.size(); t++)
936 variance += pow(log10(indVec.at(t))-obsMean,2)/numSpec;
939 for(int k = 0; k < indVec.size(); k++)
940 rO += log10(indVec.at(k));
942 double veilLine = .5;//The desired veil line.
943 double auxFunc = -(obsMean-rO)/(obsMean-log10(veilLine));
944 double uX = obsMean-auxFunc*(obsMean-log10(veilLine));
945 double vX = variance + auxFunc*pow(obsMean-log10(veilLine),2);
946 double z = (log10(veilLine)-uX)/pow(vX, .5);
947 double p = .5*(erf(z)+1);
948 double specRichness = numSpec/(1-p);
950 double numUnseen = .5*(erf((log10(.5)-uX)/pow(vX,.5))+1)*specRichness;
954 for(int i = 1; i < 8; i++)
956 double a = pow(10, i-1)+.5;
958 double c = (b - uX)/pow(vX,.5);
959 double d = .5*(erf(c)+1);
960 double numS = d*specRichness;
961 double toPush = numS - numUnseen;
962 cExp.push_back(toPush);
966 for(int i = 0; i < 8; i++)
969 for(int r = 0; r < indVec.size(); r++)
971 if(indVec.at(r) < pow(10, i-1)+.5)
972 sumOct += specVec.at(r);
975 cObs.push_back(sumOct);
976 sumOct = specVec.at(r);
979 if(r == indVec.size()-1)
980 cObs.push_back(sumOct);
984 //Statistical Analysis
985 double d = vecCalc.findDStat(cExp, cObs, numSpec);
986 cout << "DStat = " << d << "\n";
987 cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
988 cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
990 /***********************************************************************/