Changeset 3743 for trunk/yat/statistics/Fisher.cc
 Timestamp:
 Jul 12, 2018, 2:43:25 AM (4 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk
 Property svn:mergeinfo changed
/branches/0.15stable merged: 3714,3737,37393741
 Property svn:mergeinfo changed

trunk/yat/statistics/Fisher.cc
r3455 r3743 5 5 Copyright (C) 2005 Peter Johansson 6 6 Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson 7 Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015 Peter Johansson7 Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2018 Peter Johansson 8 8 9 9 This file is part of the yat library, http://dev.thep.lu.se/yat … … 87 87 88 88 /* 89 (n) n! 90 ( ) =  91 (k) k!(nk)! 92 93 94 89 95 ( n ) 90 (k+1) (k+1)!(nk)!91  =  92 (n) k!(nk1)!96 (k+1) k!(nk)! 97  =  = (nk) / (k+1) 98 (n) (k+1)!(nk1)! 93 99 (k) 94 100 */ 95 101 double Fisher::choose_ratio(unsigned int n, unsigned int k) const 96 102 { 97 if (k==0) 98 return n; 99 if (k==n1) 100 return k+1; 103 assert(k+1 <= n); 101 104 return static_cast<double>(nk) / (k+1); 102 105 } … … 106 109 unsigned int n2, unsigned int t) const 107 110 { 111 // P(X = k+1) / P(X = k) = 112 // (choose(n1, k+1) * choose(n2, tk1) / choose(n1+n2, t) ) / 113 // (choose(n1, k) * choose(n2, tk) / choose(n1+n2, t) ) = 114 // choose(n1, k+1) / choose(n1, k) /( choose(n2, tk)/choose(n2,tk1) ) = 115 // choose(n1,k+1)/choose(n1,k) / (choose(n2,(tk1)+1)/choose(n2,tk1)) 116 // choose_ratio(n1,k) / choose_ratio(n2, tk1) 117 // where choose_ratio(a, b) = choose(a, b+1) / choose(a,b) 108 118 return choose_ratio(n1, k) / choose_ratio(n2, tk1); 109 119 } … … 240 250 { 241 251 // The hypergeometric distribution is unimodal (only one peak) and the 242 // peak (mode) occurs at: (t+1)*(n1+1)/(n+1) 243 244 double mode = static_cast<double>(a_+c_+1)*static_cast<double>(a_+b_+1)/ 245 static_cast<double>(a_+b_+c_+d_+1); 252 // peak (mode) occurs at: floor((t+1)*(n1+1)/(n+1)) 253 254 double mode = std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0)); 246 255 if (a_ >= mode) 247 256 return p_value_exact(a_, a_+b_, c_+d_, a_+c_); … … 267 276 assert(n1); 268 277 assert(n2); 278 279 // special case when k=0 and mode (peak of distribution) is at 280 // zero as well. If mode is not zero 2x2 table should have been 281 // mirrored before calling this function. 282 if (k == 0) { 283 assert(std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0)) == 0.0); 284 return 1.0; 285 } 269 286 270 287 assert(k);
Note: See TracChangeset
for help on using the changeset viewer.