5 * Created by Sarah Westcott on 1/7/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "jackknife.h"
12 /***********************************************************************/
13 void Jackknife::getAMatrix(void){
15 vector<vector<double> > B = binomial(maxOrder);
17 aMat.resize(maxOrder+1);
19 for(int i=0;i<=maxOrder;i++){
20 aMat[i].resize(maxOrder+1);
21 for(int j=1;j<=maxOrder;j++){
23 aMat[i][j] = 1 + B[i][j] * (int)(pow(-1.0,j+1));
28 cout << "Standard Error: " << e.what() << " has occurred in the Jackknife class Function getAMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
32 cout << "An unknown error has occurred in the Jackknife class function getAMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
37 /**************************************************************************************************/
39 double Jackknife::CN(double z){
41 if(z>6.0) { return 0.0; }
42 if(z<-6.0) { return 0.0; }
44 const double b1= 0.31938153;
45 const double b2= -0.356563782;
46 const double b3= 1.781477937;
47 const double b4= -1.821255978;
48 const double b5= 1.330274429;
49 const double p= 0.2316419;
50 const double c2= 0.3989423;
53 double t=1.0/(1.0+a*p);
54 double b=c2*exp((-z)*(z/2.0));
55 double n=((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
60 cout << "Standard Error: " << e.what() << " has occurred in the Jackknife class Function CN. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
64 cout << "An unknown error has occurred in the Jackknife class function CN. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
69 /***********************************************************************/
71 EstOutput Jackknife::getValues(SAbundVector* rank){
73 //EstOutput jackData(3,0);
76 double jack, jacklci, jackhci;
78 double maxRank = (double)rank->getMaxRank();
79 int S = rank->getNumBins();
82 double variance[maxOrder+1];
87 for(int i=0;i<=maxOrder;i++){
90 for(int j=1;j<=maxRank;j++){
92 N[i] += aMat[i][j]*rank->get(j);
93 variance[i] += aMat[i][j]*aMat[i][j]*rank->get(j);
97 variance[i] += rank->get(j);
100 variance[i] = variance[i]-N[i];
103 for(int j=1;j<=maxRank;j++){
104 if(j<=i){ var += rank->get(j)*pow((aMat[i][j]-aMat[i-1][j]),2.0); }
105 else { var += 0.0000; }
107 var -= ((N[i]-N[i-1])*(N[i]-N[i-1]))/S;
108 var = var * S / (S-1);
109 double T = (N[i]-N[i-1])/sqrt(var);
110 if(T<=0.00){ p[i-1] = 1.00000; }
111 else{ p[i-1] = CN(T); }
118 if(i == maxOrder){ k=1; }
124 double c = (0.05-p[k-1])/(p[k]-p[k-1]);
126 jack = c*N[k]+(1-c)*N[k-1];
127 for(int j=1;j<=maxRank;j++){
128 if(j<=k){ ci += rank->get(j)*pow((c*aMat[k][j]+(1-c)*aMat[k-1][j]),2.0); }
129 else { ci += rank->get(j); }
131 ci = 1.96 * sqrt(ci - jack);
135 ci = 1.96*sqrt(variance[1]);
148 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
149 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
150 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
154 catch(exception& e) {
155 cout << "Standard Error: " << e.what() << " has occurred in the Jackknife class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
159 cout << "An unknown error has occurred in the Jackknife class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
164 /***********************************************************************/