'''See http://en.wikipedia.org/wiki/Fisher's_exact_test''' The test is used to examine the significance of the association between two variables in a 2 x 2 contingency table. # Created By : Dr. Detlef Groth # Created : Mon Mar 31 09:50:42 2008 # Last Modified : <080403.2128> # # based on code of http://www.perlmonks.org/?node_id=482919 tlm # and Alan Miller's FORTRAN-implementation http://lib.stat.cmu.edu/apstat/245 # ############################################################################## namespace eval fisher { set Tolerance 1 while {[expr {1 + $Tolerance/2.0}] > 1} { set Tolerance [expr {$Tolerance / 2.0}] } } proc fisher::test {a b c d {alternative greater}} { if {$alternative eq "less"} { return [fisherExactTest $c $d $a $b] } else { return [fisherExactTest $a $b $c $d] } } proc fisher::fisherExactTest {a b c d {ts false}} { set test [expr $a*( $a + $b + $c + $d ) - ( $a + $b )*( $a + $c )] if {$test < 0 && $ts} { return 1 } if { $test < 0 } { if { $d < $a } { set p_val [_fishers_exact $d $c $b $a $ts true] } else { set p_val [_fishers_exact $a $b $c $d $ts true] } } else { if { $b < $c } { set p_val [_fishers_exact $b $a $d $c $ts false] } else { set p_val [_fishers_exact $c $d $a $b $ts false] } } return $p_val; } proc fisher::_fishers_exact {a b c d ts complement} { foreach {aa bb cc dd} [list $a $b $c $d] { break } set delta [_single_term $aa $bb $cc $dd] set first $delta set p_val 0 while {true} { set p_val [expr {$p_val + $delta}] if {$aa < 1} { break } incr bb incr cc set delta [expr {$delta * ( 1.0* ( $aa * $dd )/( $bb * $cc ) )}] incr aa -1 incr dd -1 } if { $ts } { if {$b < $c} { set m $b } else { set m $c } set aa [expr $a+$m] set bb [expr $b-$m] set cc [expr $c-$m] set dd [expr $d+$m] set delta [_single_term $aa $bb $cc $dd] set bound [expr {-1.0*$Tolerance}] while { $bound <= [expr {( $first - $delta )/$first}] && $aa > $a } { set p_val [expr {$p_val + $delta}]; incr bb incr cc set delta [expr {$delta * ( ( $aa * $dd )/( $bb * $cc ) )}] incr aa -1 incr dd -1 } } elseif { $complement } { set p_val [expr {1 - $p_val + $first}] } return $p_val } proc fisher::_single_term {a b c d} { set r1 [expr {$a + $b}] set r2 [expr {$c + $d}] set c1 [expr {$a + $c}] set c2 [expr {$b + $d}] set N [expr {$r1 + $r2}] set s1 [expr {[_ln_fact $r1]+[_ln_fact $r2]+[_ln_fact $c1]+[_ln_fact $c2]-[_ln_fact $N]}] set s2 [expr {[_ln_fact $a]+[_ln_fact $b]+[_ln_fact $c]+[_ln_fact $d]}] return [expr {exp($s1-$s2)}] } proc fisher::_ln_fact {n} { if {$n <=1} { return 0 } return [_ln_gamm [incr n]] } proc fisher::_ln_gamm {z} { # Reference: "Lanczos, C. 'A precision approximation # of the gamma function', J. SIAM Numer. Anal., B, 1, 86-96, 1964." # Translation of Alan Miller's FORTRAN-implementation # See http://lib.stat.cmu.edu/apstat/245 set x 0 set x [expr {$x + 0.1659470187408462e-06/($z+7)}] set x [expr {$x + 0.9934937113930748e-05/($z+6)}] set x [expr {$x - 0.1385710331296526 /($z+5)}] set x [expr {$x + 12.50734324009056 /($z+4)}] set x [expr {$x - 176.6150291498386 /($z+3)}] set x [expr {$x + 771.3234287757674 /($z+2)}] set x [expr {$x - 1259.139216722289 /($z+1)}] set x [expr {$x + 676.5203681218835 /($z)}] set x [expr { $x + 0.9999999999995183}] return [expr {log($x)-5.58106146679532777-$z+($z-0.5)*log($z+6.5)}] } '''Test''' (dgroth) 127 % fisher::test 2 10 15 3 0.9999845190186861 (dgroth) 128 % fisher::test 2 10 15 3 less 0.00046518094336290726 (dgroth) 129 % fisher::test 2 10 15 3 greater 0.9999845190186861 (dgroth) 130 % fisher::test 3 1 1 3 greater 0.2428571428571418 Or the solve the first example from http://mathworld.wolfram.com/FishersExactTest.html (dgroth) 131 % fisher::test 5 1 0 4 greater 0.02380952380952382 ---- !!!!!! %| [Statistics] |% !!!!!!