Class: ViralSeq::Rubystats::FishersExactTest

Inherits:
Object
  • Object
show all
Defined in:
lib/viral_seq/rubystats.rb

Overview

Fisher’s exact test

Instance Method Summary collapse

Constructor Details

#initialize ⇒ FishersExactTest

Returns a new instance of FishersExactTest.



12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 
# File 'lib/viral_seq/rubystats.rb', line 12

def initialize
  @sn11    = 0.0
  @sn1_    = 0.0
  @sn_1    = 0.0
  @sn      = 0.0
  @sprob   = 0.0

  @sleft   = 0.0
  @sright  = 0.0
  @sless   = 0.0
  @slarg   = 0.0

  @left    = 0.0
  @right   = 0.0
  @twotail = 0.0
end

Instance Method Details

#calculate(n11_, n12_, n21_, n22_) ⇒ Object 



154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
 
# File 'lib/viral_seq/rubystats.rb', line 154

def calculate(n11_,n12_,n21_,n22_)
  n11_ *= -1 if n11_ < 0
  n12_ *= -1 if n12_ < 0
  n21_ *= -1 if n21_ < 0
  n22_ *= -1 if n22_ < 0
  n1_     = n11_ + n12_
  n_1     = n11_ + n21_
  n       = n11_ + n12_ + n21_ + n22_
  exact(n11_,n1_,n_1,n)
  left    = @sless
  right   = @slarg
  twotail = @sleft + @sright
  twotail = 1 if twotail > 1
  values_hash = { :left =>left, :right =>right, :twotail =>twotail }
  return values_hash
end

#exact(n11, n1_, n_1, n) ⇒ Object 



91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
 
# File 'lib/viral_seq/rubystats.rb', line 91

def exact(n11,n1_,n_1,n)

  p = i = j = prob = 0.0

  max = n1_
  max = n_1 if n_1 < max
  min = n1_ + n_1 - n
  min = 0 if min < 0

  if min == max
    @sless  = 1
    @sright = 1
    @sleft  = 1
    @slarg  = 1
    return 1
  end

  prob = hyper0(n11,n1_,n_1,n)
  @sleft = 0

  p = hyper(min)
  i = min + 1
  while p < (0.99999999 * prob)
    @sleft += p
    p = hyper(i)
    i += 1
  end

  i -= 1

  if p < (1.00000001*prob)
    @sleft += p
  else
    i -= 1
  end

  @sright = 0

  p = hyper(max)
  j = max - 1
  while p < (0.99999999 * prob)
    @sright += p
    p = hyper(j)
    j -= 1
  end
  j += 1

  if p < (1.00000001*prob)
    @sright += p
  else
    j += 1
  end

  if (i - n11).abs < (j - n11).abs
    @sless = @sleft
    @slarg = 1 - @sleft + prob
  else
    @sless = 1 - @sright + prob
    @slarg = @sright
  end
  return prob
end

#hyper(n11) ⇒ Object 



62
63
64
 
# File 'lib/viral_seq/rubystats.rb', line 62

def hyper(n11)
  return hyper0(n11, 0, 0, 0)
end

#hyper0(n11i, n1_i, n_1i, ni) ⇒ Object 



66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
 
# File 'lib/viral_seq/rubystats.rb', line 66

def hyper0(n11i,n1_i,n_1i,ni)
  if n1_i == 0 and n_1i ==0 and ni == 0
    unless n11i % 10 == 0
      if n11i == @sn11+1
        @sprob *= ((@sn1_ - @sn11)/(n11i.to_f))*((@sn_1 - @sn11)/(n11i.to_f + @sn - @sn1_ - @sn_1))
        @sn11 = n11i
        return @sprob
      end
      if n11i == @sn11-1
        @sprob *= ((@sn11)/(@sn1_-n11i.to_f))*((@sn11+@sn-@sn1_-@sn_1)/(@sn_1-n11i.to_f))
        @sn11 = n11i
        return @sprob
      end
    end
    @sn11 = n11i
  else
    @sn11 = n11i
    @sn1_ = n1_i
    @sn_1 = n_1i
    @sn   = ni
  end
  @sprob = hyper_323(@sn11,@sn1_,@sn_1,@sn)
  return @sprob
end

#hyper_323(n11, n1_, n_1, n) ⇒ Object 



58
59
60
 
# File 'lib/viral_seq/rubystats.rb', line 58

def hyper_323(n11, n1_, n_1, n)
  return ::Math.exp(lnbico(n1_, n11) + lnbico(n-n1_, n_1-n11) - lnbico(n, n_1))
end

#lnbico(n, k) ⇒ Object 



54
55
56
 
# File 'lib/viral_seq/rubystats.rb', line 54

def lnbico(n,k)
  return lnfact(n) - lnfact(k) - lnfact(n-k)
end

#lnfact(n) ⇒ Object 



46
47
48
49
50
51
52
 
# File 'lib/viral_seq/rubystats.rb', line 46

def lnfact(n)
  if n <= 1
    return 0
  else
    return lngamm(n+1)
  end
end

#lngamm(z) ⇒ Object

See Also:



31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
# File 'lib/viral_seq/rubystats.rb', line 31

def lngamm(z)
  x = 0
  x += 0.0000001659470187408462 / (z+7)
  x += 0.000009934937113930748  / (z+6)
  x -= 0.1385710331296526       / (z+5)
  x += 12.50734324009056        / (z+4)
  x -= 176.6150291498386        / (z+3)
  x += 771.3234287757674        / (z+2)
  x -= 1259.139216722289        / (z+1)
  x += 676.5203681218835        / (z)
  x += 0.9999999999995183

  return(::Math.log(x)-5.58106146679532777-z+(z-0.5) * ::Math.log(z+6.5))
end