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 = m->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 m->errorOut(e, "Jackknife", "getAMatrix");
33 /**************************************************************************************************/
35 double Jackknife::CN(double z){
37 if(z>6.0) { return 0.0; }
38 if(z<-6.0) { return 0.0; }
40 const double b1= 0.31938153;
41 const double b2= -0.356563782;
42 const double b3= 1.781477937;
43 const double b4= -1.821255978;
44 const double b5= 1.330274429;
45 const double p= 0.2316419;
46 const double c2= 0.3989423;
49 double t=1.0/(1.0+a*p);
50 double b=c2*exp((-z)*(z/2.0));
51 double n=((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
56 m->errorOut(e, "Jackknife", "CN");
61 /***********************************************************************/
63 EstOutput Jackknife::getValues(SAbundVector* rank){
65 //EstOutput jackData(3,0);
68 double jack, jacklci, jackhci;
70 int maxRank = (double)rank->getMaxRank();
71 int S = rank->getNumBins();
74 double variance[maxOrder+1];
79 for(int i=0;i<=maxOrder;i++){
82 for(int j=1;j<=maxRank;j++){
84 N[i] += aMat[i][j]*rank->get(j);
85 variance[i] += aMat[i][j]*aMat[i][j]*rank->get(j);
89 variance[i] += rank->get(j);
92 variance[i] = variance[i]-N[i];
95 for(int j=1;j<=maxRank;j++){
96 if(j<=i){ var += rank->get(j)*pow((aMat[i][j]-aMat[i-1][j]),2.0); }
97 else { var += 0.0000; }
99 var -= ((N[i]-N[i-1])*(N[i]-N[i-1]))/S;
100 var = var * S / (S-1);
101 double T = (N[i]-N[i-1])/sqrt(var);
102 if(T<=0.00){ p[i-1] = 1.00000; }
103 else{ p[i-1] = CN(T); }
110 if(i == maxOrder){ k=1; }
116 double c = (0.05-p[k-1])/(p[k]-p[k-1]);
118 jack = c*N[k]+(1-c)*N[k-1];
119 for(int j=1;j<=maxRank;j++){
120 if(j<=k){ ci += rank->get(j)*pow((c*aMat[k][j]+(1-c)*aMat[k-1][j]),2.0); }
121 else { ci += rank->get(j); }
123 ci = 1.96 * sqrt(ci - jack);
127 ci = 1.96*sqrt(variance[1]);
140 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
141 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
142 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
146 catch(exception& e) {
147 m->errorOut(e, "Jackknife", "getValues");
152 /***********************************************************************/