5 * Created by westcott on 7/8/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 //translated to c++ using source code http://www.langsrud.com/stat/fisher.htm as a reference
12 #include "mothurfisher.h"
13 /***********************************************************/
14 double MothurFisher::fexact(double n11_, double n12_, double n21_, double n22_) {
16 sleft = 0.0; sright = 0.0; sless = 0.0; slarg = 0.0;
18 if(n11_<0) n11_ *= -1;
19 if(n12_<0) n12_ *= -1;
20 if(n21_<0) n21_ *= -1;
21 if(n22_<0) n22_ *= -1;
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;
29 if(twotail>1) twotail=1;
30 double result = twotail;
33 }catch(exception& e) {
34 m->errorOut(e, "MothurFisher", "fexact");
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
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;
56 return(log(x)-5.58106146679532777-z+(z-0.5)*log(z+6.5));
58 }catch(exception& e) {
59 m->errorOut(e, "MothurFisher", "lngamm");
64 /***********************************************************/
65 double MothurFisher::lnfact(double n){
69 }catch(exception& e) {
70 m->errorOut(e, "MothurFisher", "lnfact");
74 /***********************************************************/
75 double MothurFisher::lnbico(double n, double k){
77 return(lnfact(n)-lnfact(k)-lnfact(n-k));
78 }catch(exception& e) {
79 m->errorOut(e, "MothurFisher", "lnbico");
83 /***********************************************************/
84 double MothurFisher::hyper_323(double n11, double n1_, double n_1, double n){
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");
92 /***********************************************************/
93 double MothurFisher::myhyper(double n11){
95 double hyper0Result = hyper0(n11,0,0,0);
97 }catch(exception& e) {
98 m->errorOut(e, "MothurFisher", "myhyper");
102 /***********************************************************/
103 double MothurFisher::hyper0(double n11i, double n1_i, double n_1i, double ni) {
105 if (!((n1_i != 0)&&(n_1i != 0)&&(ni != 0))) {
106 if(!(((int)n11i % 10) == 0)){
109 sprob *= ((sn1_-sn11)/(n11i))*((sn_1-sn11)/(n11i+sn-sn1_-sn_1));
115 sprob *= ((sn11)/(sn1_-n11i))*((sn11+sn-sn1_-sn_1)/(sn_1-n11i));
128 sprob = hyper_323(sn11,sn1_,sn_1,sn);
131 }catch(exception& e) {
132 m->errorOut(e, "MothurFisher", "hyper0");
136 /***********************************************************/
137 double MothurFisher::exact(double n11, double n1_, double n_1, double n){
142 double min = n1_+n_1-n;
152 prob=hyper0(n11,n1_,n_1,n);
155 for(i=min+1; p<0.99999999*prob; i++)
161 if(p<1.00000001*prob) sleft += p;
165 for(j=max-1; p<0.99999999*prob; j--)
171 if(p<1.00000001*prob) sright += p;
173 if(abs(i-n11)<abs(j-n11))
176 slarg = 1 - sleft + prob;
180 sless = 1 - sright + prob;
185 }catch(exception& e) {
186 m->errorOut(e, "MothurFisher", "exact");
190 /***********************************************************/