cistematic/cisstat/fisher.py

   1 ###########################################################################
   2 #                                                                         #
   3 # C O P Y R I G H T   N O T I C E                                         #
   4 #  Copyright (c) 2003-10 by:                                              #
   5 #    * California Institute of Technology                                 #
   6 #                                                                         #
   7 #    All Rights Reserved.                                                 #
   8 #                                                                         #
   9 # Permission is hereby granted, free of charge, to any person             #
  10 # obtaining a copy of this software and associated documentation files    #
  11 # (the "Software"), to deal in the Software without restriction,          #
  12 # including without limitation the rights to use, copy, modify, merge,    #
  13 # publish, distribute, sublicense, and/or sell copies of the Software,    #
  14 # and to permit persons to whom the Software is furnished to do so,       #
  15 # subject to the following conditions:                                    #
  16 #                                                                         #
  17 # The above copyright notice and this permission notice shall be          #
  18 # included in all copies or substantial portions of the Software.         #
  19 #                                                                         #
  20 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
  21 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
  22 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
  23 # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
  24 # BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
  25 # ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
  26 # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
  27 # SOFTWARE.                                                               #
  28 ###########################################################################
  29 #
  30 #
  31 # Fisher's exact test, as described in Gopal Kanji's 100 Statistical Tests
  32 # Table has form:
  33 #                               A   |    B      | A + B
  34 #                           -----------------------------
  35 #                               C   |    D      | C + D
  36 #                           -----------------------------
  37 #                           A + C | B + D   | A + B + C + D
  38 #
  39
  40 from math import exp
  41 from helper import factln
  42
  43 def fisher(a, b, c, d):
  44     lowest = 0
  45     if a < b and a < c and a < d:
  46         lowest = "a"
  47         lowestNum = a
  48     elif b < c and b < d:
  49         lowest = "b"
  50         lowestNum = b
  51     elif c < d:
  52         lowest = "c"
  53         lowestNum = c
  54     else:
  55         lowest = "d"
  56         lowestNum = d
  57
  58     origAB = a + b
  59     origCD = c + d
  60     origAC = a + c
  61     origBD = b + d
  62     N = a + b + c + d
  63
  64     first = factln(origAB) + factln(origCD) + factln(origAC) + factln(origBD) - factln(N)
  65
  66     prob = 0.0
  67     for x in range(lowestNum + 1):
  68         if lowest == "a":
  69             newA = x
  70             newC = origAC - newA
  71             newB = origAB - newA
  72             newD = origCD - newC
  73         elif lowest == "b":
  74             newB = x
  75             newA = origAB - newB
  76             newD = origBD - newB
  77             newC = origAC - newA
  78         elif lowest == "c":
  79             newC = x
  80             newA = origAC - newC
  81             newD = origCD - newC
  82             newB = origAB - newA
  83         else:
  84             newD = x
  85             newB = origBD - newD
  86             newC = origCD - newD
  87             newA = origAB - newB
  88
  89         currentProb =  exp(first - ((factln(newA) + factln(newB) + factln(newC) + factln(newD))))
  90         #print 'did %f %f %f %f: %f' % (newA, newB, newC, newD, currentProb)
  91         prob += currentProb
  92
  93     if prob > 1.0:
  94         prob = 1.0
  95
  96     return prob