mothurfisher.cpp

   1 /*
   2  *  mothurfisher.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 7/8/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 //translated to c++ using source code http://www.langsrud.com/stat/fisher.htm as a reference
  11
  12 #include "mothurfisher.h"
  13 /***********************************************************/
  14 double MothurFisher::fexact(double n11_, double n12_, double n21_, double n22_) {
  15         try {
  16                 sleft = 0.0; sright = 0.0; sless = 0.0; slarg = 0.0;
  17
  18                 if(n11_<0) n11_ *= -1;
  19                 if(n12_<0) n12_ *= -1;
  20                 if(n21_<0) n21_ *= -1;
  21                 if(n22_<0) n22_ *= -1;
  22
  23                 double n1_ = n11_+n12_;
  24                 double n_1 = n11_+n21_;
  25                 double n   = n11_ +n12_ +n21_ +n22_;
  26                 exact(n11_,n1_,n_1,n);
  27                 double twotail = sleft+sright;
  28
  29                 if(twotail>1) twotail=1;
  30                 double result = twotail;
  31                 return result;
  32
  33         }catch(exception& e) {
  34                 m->errorOut(e, "MothurFisher", "fexact");
  35                 exit(1);
  36         }
  37 }
  38 /***********************************************************/
  39 double MothurFisher::lngamm(double z) {
  40         // Reference: "Lanczos, C. 'A precision approximation
  41         // of the gamma function', J. SIAM Numer. Anal., B, 1, 86-96, 1964."
  42         // Translation of  Alan Miller's FORTRAN-implementation
  43         // See http://lib.stat.cmu.edu/apstat/245
  44         try {
  45                 double x = 0;
  46                 x += 0.1659470187408462e-06/(z+7);
  47                 x += 0.9934937113930748e-05/(z+6);
  48                 x -= 0.1385710331296526    /(z+5);
  49                 x += 12.50734324009056     /(z+4);
  50                 x -= 176.6150291498386     /(z+3);
  51                 x += 771.3234287757674     /(z+2);
  52                 x -= 1259.139216722289     /(z+1);
  53                 x += 676.5203681218835     /(z);
  54                 x += 0.9999999999995183;
  55
  56                 return(log(x)-5.58106146679532777-z+(z-0.5)*log(z+6.5));
  57
  58         }catch(exception& e) {
  59                 m->errorOut(e, "MothurFisher", "lngamm");
  60                 exit(1);
  61         }
  62 }
  63
  64 /***********************************************************/
  65 double MothurFisher::lnfact(double n){
  66         try {
  67                 if(n <= 1) return(0);
  68                 return(lngamm(n+1));
  69         }catch(exception& e) {
  70                 m->errorOut(e, "MothurFisher", "lnfact");
  71                 exit(1);
  72         }
  73 }
  74 /***********************************************************/
  75 double MothurFisher::lnbico(double n, double k){
  76         try {
  77                 return(lnfact(n)-lnfact(k)-lnfact(n-k));
  78         }catch(exception& e) {
  79                 m->errorOut(e, "MothurFisher", "lnbico");
  80                 exit(1);
  81         }
  82 }
  83 /***********************************************************/
  84 double MothurFisher::hyper_323(double n11, double n1_, double n_1, double n){
  85         try {
  86                 return(exp(lnbico(n1_,n11)+lnbico(n-n1_,n_1-n11)-lnbico(n,n_1)));
  87         }catch(exception& e) {
  88                 m->errorOut(e, "MothurFisher", "hyper_323");
  89                 exit(1);
  90         }
  91 }
  92 /***********************************************************/
  93 double MothurFisher::myhyper(double n11){
  94         try {
  95                 double hyper0Result = hyper0(n11,0,0,0);
  96                 return hyper0Result;
  97         }catch(exception& e) {
  98                 m->errorOut(e, "MothurFisher", "myhyper");
  99                 exit(1);
 100         }
 101 }
 102 /***********************************************************/
 103 double MothurFisher::hyper0(double n11i, double n1_i, double n_1i, double ni) {
 104         try {
 105                 if (!((n1_i != 0)&&(n_1i != 0)&&(ni != 0))) {
 106                         if(!(((int)n11i % 10) == 0)){
 107                                 if(n11i==sn11+1)
 108                                 {
 109                                         sprob *= ((sn1_-sn11)/(n11i))*((sn_1-sn11)/(n11i+sn-sn1_-sn_1));
 110                                         sn11 = n11i;
 111                                         return sprob;
 112                                 }
 113                                 if(n11i==sn11-1)
 114                                 {
 115                                         sprob *= ((sn11)/(sn1_-n11i))*((sn11+sn-sn1_-sn_1)/(sn_1-n11i));
 116                                         sn11 = n11i;
 117                                         return sprob;
 118                                 }
 119                         }
 120                         sn11 = n11i;
 121                 }else{
 122                         sn11 = n11i;
 123                         sn1_=n1_i;
 124                         sn_1=n_1i;
 125                         sn=ni;
 126                 }
 127
 128                 sprob = hyper_323(sn11,sn1_,sn_1,sn);
 129                 return sprob;
 130
 131         }catch(exception& e) {
 132                 m->errorOut(e, "MothurFisher", "hyper0");
 133                 exit(1);
 134         }
 135 }
 136 /***********************************************************/
 137 double MothurFisher::exact(double n11, double n1_, double n_1, double n){
 138         try {
 139                 double p,i,j,prob;
 140                 double max=n1_;
 141                 if(n_1<max) max=n_1;
 142                 double min = n1_+n_1-n;
 143                 if(min<0) min=0;
 144                 if(min==max)
 145                 {
 146                         sless = 1;
 147                         sright= 1;
 148                         sleft = 1;
 149                         slarg = 1;
 150                         return 1;
 151                 }
 152                 prob=hyper0(n11,n1_,n_1,n);
 153                 sleft=0;
 154                 p=myhyper(min);
 155                 for(i=min+1; p<0.99999999*prob; i++)
 156                 {
 157                         sleft += p;
 158                         p=myhyper(i);
 159                 }
 160                 i--;
 161                 if(p<1.00000001*prob) sleft += p;
 162                 else i--;
 163                 sright=0;
 164                 p=myhyper(max);
 165                 for(j=max-1; p<0.99999999*prob; j--)
 166                 {
 167                         sright += p;
 168                         p=myhyper(j);
 169                 }
 170                 j++;
 171                 if(p<1.00000001*prob) sright += p;
 172                 else j++;
 173                 if(abs(i-n11)<abs(j-n11))
 174                 {
 175                         sless = sleft;
 176                         slarg = 1 - sleft + prob;
 177                 }
 178                 else
 179                 {
 180                         sless = 1 - sright + prob;
 181                         slarg = sright;
 182                 }
 183                 return prob;
 184
 185         }catch(exception& e) {
 186                 m->errorOut(e, "MothurFisher", "exact");
 187                 exit(1);
 188         }
 189 }
 190 /***********************************************************/
 191
 192
 193