## FisherTest

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
#
# 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

 Category Statistics Category Mathematics