test/test_Prob.py

   1 #!/usr/bin/env python
   2 # Time-stamp: <2019-09-25 11:35:54 taoliu>
   3
   4 """Module Description: Test functions to calculate probabilities.
   5
   6 This code is free software; you can redistribute it and/or modify it
   7 under the terms of the BSD License (see the file LICENSE included with
   8 the distribution).
   9 """
  10
  11 import unittest
  12
  13 from math import log10
  14 from MACS2.Prob import *
  15
  16 # ------------------------------------
  17 # Main function
  18 # ------------------------------------
  19
  20 class Test_factorial(unittest.TestCase):
  21
  22     def setUp(self):
  23         self.n1 = 100
  24         self.n2 = 10
  25         self.n3 = 1
  26
  27     def test_factorial_big_n1(self):
  28         expect = 9.332622e+157
  29         result = factorial(self.n1)
  30         self.assertTrue( abs(result - expect) < 1e-5*result)
  31
  32     def test_factorial_median_n2(self):
  33         expect = 3628800
  34         result = factorial(self.n2)
  35         self.assertEqual(result, expect)
  36
  37     def test_factorial_small_n3(self):
  38         expect = 1
  39         result = factorial(self.n3)
  40         self.assertEqual(result, expect)
  41
  42 class Test_poisson_cdf(unittest.TestCase):
  43
  44     def setUp(self):
  45         # n, lam
  46         self.n1 = (80,100)
  47         self.n2 = (200,100)
  48         self.n3 = (100,1000)
  49         self.n4 = (1500,1000)
  50
  51     def test_poisson_cdf_n1(self):
  52         expect = (round(0.9773508,5),round(0.02264918,5))
  53         result = (round(poisson_cdf(self.n1[0],self.n1[1],False),5),
  54                   round(poisson_cdf(self.n1[0],self.n1[1],True),5))
  55         self.assertEqual( result, expect )
  56
  57     def test_poisson_cdf_n2(self):
  58         expect = (round(log10(4.626179e-19),4),
  59                   round(log10(1),4))
  60         result = (round(log10(poisson_cdf(self.n2[0],self.n2[1],False)),4),
  61                   round(log10(poisson_cdf(self.n2[0],self.n2[1],True)),4))
  62         self.assertEqual( result, expect )
  63
  64     def test_poisson_cdf_n3(self):
  65         expect = (round(log10(1),2),
  66                   round(log10(6.042525e-293),2))
  67         result = (round(poisson_cdf(self.n3[0],self.n3[1],False,True),2),
  68                   round(poisson_cdf(self.n3[0],self.n3[1],True,True),2))
  69         self.assertEqual( result, expect )
  70
  71     def test_poisson_cdf_n4(self):
  72         expect = (round(log10(2.097225e-49),4),
  73                   round(log10(1),4))
  74         result = (round(log10(poisson_cdf(self.n4[0],self.n4[1],False)),4),
  75                   round(log10(poisson_cdf(self.n4[0],self.n4[1],True)),4))
  76         self.assertEqual( result, expect )
  77
  78 class Test_chisq_p_e(unittest.TestCase):
  79     """Test chisq pvalue calculation -- assuming df is an even number. We
  80     only implemented even number pchisq for upper tail. Because this
  81     is the function we need to combine p-values using fisher's method
  82
  83     """
  84     def setUp(self):
  85         # x, k, p(upper), -log p upper, -log10 p upper
  86         self.c = ((10, 2, 0.006737947, 5, 2.171472),
  87                   (100, 2, 1.92875e-22, 50, 21.71472),
  88                   (1000, 22, 1.956374e-197, 452.9382, 196.7085),
  89                   (10, 4, 0.04042768, 3.208241, 1.393321),
  90                   (100, 8, 4.269159e-18, 39.99511, 17.36966),
  91                   (1000, 80, 6.889598e-159, 364.181, 158.1618),
  92                   (54, 6, 7.377151e-10, 21.02746, 9.132111),
  93                   (565, 10, 5.518772e-115, 263.0891, 114.2582 ),
  94                   (7765, 12, 0, 3845.965, 1670.2814),
  95                  )
  96
  97     def test_chisq_p(self):
  98         expect = [round(x[2],4) for x in self.c]
  99         result = [round(chisq_pvalue_e(x[0],x[1]),4) for x in self.c]
 100         self.assertEqual( result, expect )
 101
 102     def test_chisq_logp(self):
 103         expect = [round(x[3],4) for x in self.c]
 104         result = [round(chisq_logp_e(x[0],x[1]),4) for x in self.c]
 105         self.assertEqual( result, expect )
 106
 107     def test_chisq_log10p(self):
 108         expect = [round(x[4],4) for x in self.c]
 109         result = [round(chisq_logp_e(x[0],x[1],log10=True),4) for x in self.c]
 110         self.assertEqual( result, expect )
 111
 112 class Test_binomial_cdf(unittest.TestCase):
 113
 114     def setUp(self):
 115         # x, a, b
 116         self.n1 = (20,1000,0.01)
 117         self.n2 = (200,1000,0.01)
 118
 119     def test_binomial_cdf_n1(self):
 120         expect = (round(0.001496482,5),round(0.9985035,5))
 121         result = (round(binomial_cdf(self.n1[0],self.n1[1],self.n1[2],False),5),
 122                   round(binomial_cdf(self.n1[0],self.n1[1],self.n1[2],True),5))
 123         self.assertEqual( result, expect )
 124
 125     def test_binomial_cdf_n2(self):
 126         expect = (round(log10(8.928717e-190),4),
 127                   round(log10(1),4))
 128         result = (round(log10(binomial_cdf(self.n2[0],self.n2[1],self.n2[2],False)),4),
 129                   round(log10(binomial_cdf(self.n2[0],self.n2[1],self.n2[2],True)),4))
 130         self.assertEqual( result, expect )
 131
 132 class Test_binomial_cdf_inv(unittest.TestCase):
 133
 134     def setUp(self):
 135         # x, a, b
 136         self.n1 = (0.1,1000,0.01)
 137         self.n2 = (0.01,1000,0.01)
 138
 139     def test_binomial_cdf_inv_n1(self):
 140         expect = 6
 141         result = binomial_cdf_inv(self.n1[0],self.n1[1],self.n1[2])
 142         self.assertEqual( result, expect )
 143
 144     def test_poisson_cdf_inv_n2(self):
 145         expect = 3
 146         result = binomial_cdf_inv(self.n2[0],self.n2[1],self.n2[2])
 147         self.assertEqual( result, expect )
 148
 149 class Test_binomial_pdf(unittest.TestCase):
 150
 151     def setUp(self):
 152         # x, a, b
 153         self.n1 = (20,1000,0.01)
 154         self.n2 = (200,1000,0.01)
 155
 156     def test_binomial_cdf_inv_n1(self):
 157         expect = round(0.001791878,5)
 158         result = round(binomial_pdf(self.n1[0],self.n1[1],self.n1[2]),5)
 159         self.assertEqual( result, expect )
 160
 161     def test_poisson_cdf_inv_n2(self):
 162         expect = round(log10(2.132196e-188),4)
 163         result = binomial_pdf(self.n2[0],self.n2[1],self.n2[2])
 164         result = round(log10(result),4)
 165         self.assertEqual( result, expect )
 166