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 = m->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                 m->errorOut(e, "Jackknife", "getAMatrix");
  29                 exit(1);
  30         }
  31 }
  32
  33 /**************************************************************************************************/
  34
  35 double Jackknife::CN(double z){
  36         try {
  37                 if(z>6.0)       {       return 0.0;             }
  38                 if(z<-6.0)      {       return 0.0;             }
  39
  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;
  47
  48                 double a=abs(z);
  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;
  52                 n = 2*b*n;
  53                 return n;
  54         }
  55         catch(exception& e) {
  56                 m->errorOut(e, "Jackknife", "CN");
  57                 exit(1);
  58         }
  59 }
  60
  61 /***********************************************************************/
  62
  63 EstOutput Jackknife::getValues(SAbundVector* rank){
  64         try {
  65                 //EstOutput jackData(3,0);
  66                 data.resize(3,0);
  67
  68                 double jack, jacklci, jackhci;
  69
  70                 int maxRank = (double)rank->getMaxRank();
  71                 int S = rank->getNumBins();
  72
  73                 double N[maxOrder+1];
  74                 double variance[maxOrder+1];
  75                 double p[maxOrder+1];
  76
  77                 int k = 0;
  78
  79                 for(int i=0;i<=maxOrder;i++){
  80                         N[i]=0.0000;
  81                         variance[i]=0.0000;
  82                         for(int j=1;j<=maxRank;j++){
  83                                 if(j<=i){
  84                                         N[i] += aMat[i][j]*rank->get(j);
  85                                         variance[i] += aMat[i][j]*aMat[i][j]*rank->get(j);
  86                                 }
  87                                 else{
  88                                         N[i] += rank->get(j);
  89                                         variance[i] += rank->get(j);
  90                                 }
  91                         }
  92                         variance[i] = variance[i]-N[i];
  93                         double var = 0.0000;
  94                         if(i>0){
  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;  }
  98                                 }
  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);                 }
 104
 105                                 if(p[i-1]>=0.05){
 106                                         k = i-1;
 107                                         break;
 108                                 }
 109                         }
 110                         if(i == maxOrder){      k=1;    }
 111                 }
 112
 113                 double ci = 0;
 114
 115                 if(k>1){
 116                         double c = (0.05-p[k-1])/(p[k]-p[k-1]);
 117                         ci = 0.0000;
 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);     }
 122                         }
 123                         ci = 1.96 * sqrt(ci - jack);
 124                 }
 125                 else if(k==1){
 126                         jack = N[1];
 127                         ci = 1.96*sqrt(variance[1]);
 128                 }else{
 129                         jack = 0.0;
 130                         ci = 0.0;
 131                 }
 132
 133                 jacklci = jack-ci;
 134                 jackhci = jack+ci;
 135
 136                 data[0] = jack;
 137                 data[1] = jacklci;
 138                 data[2] = jackhci;
 139
 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; }
 143
 144                 return data;
 145         }
 146         catch(exception& e) {
 147                 m->errorOut(e, "Jackknife", "getValues");
 148                 exit(1);
 149         }
 150 }
 151
 152 /***********************************************************************/