Class: Rubystats::FishersExactTest

Inherits:

Object

• Object
• Rubystats::FishersExactTest

Defined in:

lib/viral_seq/math.rb

Instance Method Summary

• #calculate(n11_, n12_, n21_, n22_) ⇒ Object
• #exact(n11, n1_, n_1, n) ⇒ Object
• #hyper(n11) ⇒ Object
• #hyper0(n11i, n1_i, n_1i, ni) ⇒ Object
• #hyper_323(n11, n1_, n_1, n) ⇒ Object
• #initialize ⇒ FishersExactTest constructor

  A new instance of FishersExactTest.

• #lnbico(n, k) ⇒ Object
• #lnfact(n) ⇒ Object
• #lngamm(z) ⇒ Object

  Reference: “Lanczos, C.

Constructor Details

#initialize ⇒ FishersExactTest

Returns a new instance of FishersExactTest.

# File 'lib/viral_seq/math.rb', line 202

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

# File 'lib/viral_seq/math.rb', line 346

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

# File 'lib/viral_seq/math.rb', line 283

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

# File 'lib/viral_seq/math.rb', line 254

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

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

# File 'lib/viral_seq/math.rb', line 258

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

# File 'lib/viral_seq/math.rb', line 250

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

# File 'lib/viral_seq/math.rb', line 246

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

#lnfact(n) ⇒ Object

# File 'lib/viral_seq/math.rb', line 238

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

#lngamm(z) ⇒ Object

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 lib.stat.cmu.edu/apstat/245

# File 'lib/viral_seq/math.rb', line 223

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