jackknife.cpp

   1 /*
   2  *  jacknife.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/7/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "jackknife.h"
  11
  12 /***********************************************************************/
  13 void Jackknife::getAMatrix(void){
  14         try {
  15                 vector<vector<double> > B = binomial(maxOrder);
  16
  17                 aMat.resize(maxOrder+1);
  18
  19                 for(int i=0;i<=maxOrder;i++){
  20                         aMat[i].resize(maxOrder+1);
  21                         for(int j=1;j<=maxOrder;j++){
  22
  23                                 aMat[i][j] = 1 + B[i][j] * (int)(pow(-1.0,j+1));
  24                         }
  25                 }
  26         }
  27         catch(exception& e) {
  28                 cout << "Standard Error: " << e.what() << " has occurred in the Jackknife class Function getAMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  29                 exit(1);
  30         }
  31         catch(...) {
  32                 cout << "An unknown error has occurred in the Jackknife class function getAMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  33                 exit(1);
  34         }
  35 }
  36
  37 /**************************************************************************************************/
  38
  39 double Jackknife::CN(double z){
  40         try {
  41                 if(z>6.0)       {       return 0.0;             }
  42                 if(z<-6.0)      {       return 0.0;             }
  43
  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;
  51
  52                 double a=abs(z);
  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;
  56                 n = 2*b*n;
  57                 return n;
  58         }
  59         catch(exception& e) {
  60                 cout << "Standard Error: " << e.what() << " has occurred in the Jackknife class Function CN. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  61                 exit(1);
  62         }
  63         catch(...) {
  64                 cout << "An unknown error has occurred in the Jackknife class function CN. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  65                 exit(1);
  66         }
  67 }
  68
  69 /***********************************************************************/
  70
  71 EstOutput Jackknife::getValues(SAbundVector* rank){
  72         try {
  73                 //EstOutput jackData(3,0);
  74                 data.resize(3,0);
  75
  76                 double jack, jacklci, jackhci;
  77
  78                 double maxRank = (double)rank->getMaxRank();
  79                 int S = rank->getNumBins();
  80
  81                 double N[maxOrder+1];
  82                 double variance[maxOrder+1];
  83                 double p[maxOrder+1];
  84
  85                 int k = 0;
  86
  87                 for(int i=0;i<=maxOrder;i++){
  88                         N[i]=0.0000;
  89                         variance[i]=0.0000;
  90                         for(int j=1;j<=maxRank;j++){
  91                                 if(j<=i){
  92                                         N[i] += aMat[i][j]*rank->get(j);
  93                                         variance[i] += aMat[i][j]*aMat[i][j]*rank->get(j);
  94                                 }
  95                                 else{
  96                                         N[i] += rank->get(j);
  97                                         variance[i] += rank->get(j);
  98                                 }
  99                         }
 100                         variance[i] = variance[i]-N[i];
 101                         double var = 0.0000;
 102                         if(i>0){
 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;  }
 106                                 }
 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);                 }
 112
 113                                 if(p[i-1]>=0.05){
 114                                         k = i-1;
 115                                         break;
 116                                 }
 117                         }
 118                         if(i == maxOrder){      k=1;    }
 119                 }
 120
 121                 double ci = 0;
 122
 123                 if(k>1){
 124                         double c = (0.05-p[k-1])/(p[k]-p[k-1]);
 125                         ci = 0.0000;
 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);     }
 130                         }
 131                         ci = 1.96 * sqrt(ci - jack);
 132                 }
 133                 else if(k=1){
 134                         jack = N[1];
 135                         ci = 1.96*sqrt(variance[1]);
 136                 }else{
 137                         jack = 0.0;
 138                         ci = 0.0;
 139                 }
 140
 141                 jacklci = jack-ci;
 142                 jackhci = jack+ci;
 143
 144                 data[0] = jack;
 145                 data[1] = jacklci;
 146                 data[2] = jackhci;
 147
 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; }
 151
 152                 return data;
 153         }
 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";
 156                 exit(1);
 157         }
 158         catch(...) {
 159                 cout << "An unknown error has occurred in the Jackknife class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 160                 exit(1);
 161         }
 162 }
 163
 164 /***********************************************************************/