#!/usr/bin/env python """Tests of statistical probability distribution integrals. Currently using tests against calculations in R, spreadsheets being unreliable. """ from cogent.util.unit_test import TestCase, main from cogent.maths.stats.distribution import z_low, z_high, zprob, chi_low, \ chi_high, t_low, t_high, tprob, poisson_high, poisson_low, poisson_exact, \ binomial_high, binomial_low, binomial_exact, f_low, f_high, fprob, \ stdtr, bdtr, bdtrc, pdtr, pdtrc, fdtr, fdtrc, gdtr, gdtrc, chdtri, stdtri,\ pdtri, bdtri, fdtri, gdtri __author__ = "Rob Knight" __copyright__ = "Copyright 2007-2012, The Cogent Project" __credits__ = ["Gavin Huttley", "Rob Knight", "Sandra Smit"] __license__ = "GPL" __version__ = "1.5.3" __maintainer__ = "Rob Knight" __email__ = "rob@spot.colorado.edu" __status__ = "Production" class DistributionsTests(TestCase): """Tests of particular statistical distributions.""" def setUp(self): self.values = [0, 0.01, 0.1, 0.5, 1, 2, 5, 10, 20, 30, 50, 200] self.negvalues = [-i for i in self.values] self.df = [1, 10, 100] def test_z_low(self): """z_low should match R's pnorm() function""" probs = [ 0.5000000, 0.5039894, 0.5398278, 0.6914625, 0.8413447, 0.9772499, 0.9999997, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, ] negprobs = [ 5.000000e-01, 4.960106e-01, 4.601722e-01, 3.085375e-01, 1.586553e-01, 2.275013e-02, 2.866516e-07, 7.619853e-24, 2.753624e-89, 4.906714e-198, 0.000000e+00, 0.000000e+00] for z, p in zip(self.values, probs): self.assertFloatEqual(z_low(z), p) for z, p in zip(self.negvalues, negprobs): self.assertFloatEqual(z_low(z), p) def test_z_high(self): """z_high should match R's pnorm(lower.tail=FALSE) function""" negprobs = [ 0.5000000, 0.5039894, 0.5398278, 0.6914625, 0.8413447, 0.9772499, 0.9999997, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, ] probs = [ 5.000000e-01, 4.960106e-01, 4.601722e-01, 3.085375e-01, 1.586553e-01, 2.275013e-02, 2.866516e-07, 7.619853e-24, 2.753624e-89, 4.906714e-198, 0.000000e+00, 0.000000e+00] for z, p in zip(self.values, probs): self.assertFloatEqual(z_high(z), p) for z, p in zip(self.negvalues, negprobs): self.assertFloatEqual(z_high(z), p) def test_zprob(self): """zprob should match twice the z_high probability for abs(z)""" probs = [2*i for i in [ 5.000000e-01, 4.960106e-01, 4.601722e-01, 3.085375e-01, 1.586553e-01, 2.275013e-02, 2.866516e-07, 7.619853e-24, 2.753624e-89, 4.906714e-198, 0.000000e+00, 0.000000e+00]] for z, p in zip(self.values, probs): self.assertFloatEqual(zprob(z), p) for z, p in zip(self.negvalues, probs): self.assertFloatEqual(zprob(z), p) def test_chi_low(self): """chi_low should match R's pchisq() function""" probs = { 1: [ 0.00000000, 0.07965567, 0.24817037, 0.52049988, 0.68268949, 0.84270079, 0.97465268, 0.99843460, 0.99999226, 0.99999996, 1.00000000, 1.00000000 ], 10: [ 0.000000e+00, 2.593339e-14, 2.497951e-09, 6.611711e-06, 1.721156e-04, 3.659847e-03, 1.088220e-01, 5.595067e-01, 9.707473e-01, 9.991434e-01, 9.999997e-01, 1.000000e+00, ], 100: [ 0.000000e+00, 2.906006e-180, 2.780588e-130, 2.029952e-95, 1.788777e-80, 1.233751e-65, 2.238699e-46, 2.181059e-32, 1.854727e-19, 9.056126e-13, 6.953305e-06, 1.000000e-00 ], } for df in self.df: for x, p in zip(self.values, probs[df]): self.assertFloatEqual(chi_low(x, df), p) def test_chi_high(self): """chi_high should match R's pchisq(lower.tail=FALSE) function""" probs = { 1: [ 1.000000e+00, 9.203443e-01, 7.518296e-01, 4.795001e-01, 3.173105e-01, 1.572992e-01, 2.534732e-02, 1.565402e-03, 7.744216e-06, 4.320463e-08, 1.537460e-12, 2.088488e-45, ], 10: [ 1.000000e+00, 1.000000e-00, 1.000000e-00, 9.999934e-01, 9.998279e-01, 9.963402e-01, 8.911780e-01, 4.404933e-01, 2.925269e-02, 8.566412e-04, 2.669083e-07, 1.613931e-37, ], 100:[ 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 9.99993e-01, 1.17845e-08, ], } for df in self.df: for x, p in zip(self.values, probs[df]): self.assertFloatEqual(chi_high(x, df), p) def test_t_low(self): """t_low should match R's pt() function""" probs = { 1: [ 0.5000000, 0.5031830, 0.5317255, 0.6475836, 0.7500000, 0.8524164, 0.9371670, 0.9682745, 0.9840977, 0.9893936, 0.9936347, 0.9984085, ], 10: [ 0.5000000, 0.5038910, 0.5388396, 0.6860532, 0.8295534, 0.9633060, 0.9997313, 0.9999992, 1.0000000, 1.0000000, 1.0000000, 1.0000000, ], 100:[ 0.5000000, 0.5039794, 0.5397277, 0.6909132, 0.8401379, 0.9758939, 0.9999988, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, ], } negprobs = { 1: [ 0.500000000, 0.496817007, 0.468274483, 0.352416382, 0.250000000, 0.147583618, 0.062832958, 0.031725517, 0.015902251, 0.010606402, 0.006365349, 0.001591536, ], 10: [ 5.000000e-01, 4.961090e-01, 4.611604e-01, 3.139468e-01, 1.704466e-01, 3.669402e-02, 2.686668e-04, 7.947766e-07, 1.073031e-09, 1.980896e-11, 1.237155e-13, 1.200254e-19, ], 100:[ 5.000000e-01, 4.960206e-01, 4.602723e-01, 3.090868e-01, 1.598621e-01, 2.410609e-02, 1.225087e-06, 4.950844e-17, 4.997134e-37, 4.190166e-52, 7.236082e-73, 2.774197e-132, ], } for df in self.df: for x, p in zip(self.values, probs[df]): self.assertFloatEqualRel(t_low(x, df), p, eps=1e-4) for x, p in zip(self.negvalues, negprobs[df]): self.assertFloatEqualRel(t_low(x, df), p, eps=1e-4) def test_t_high(self): """t_high should match R's pt(lower.tail=FALSE) function""" negprobs = { 1: [ 0.5000000, 0.5031830, 0.5317255, 0.6475836, 0.7500000, 0.8524164, 0.9371670, 0.9682745, 0.9840977, 0.9893936, 0.9936347, 0.9984085, ], 10: [ 0.5000000, 0.5038910, 0.5388396, 0.6860532, 0.8295534, 0.9633060, 0.9997313, 0.9999992, 1.0000000, 1.0000000, 1.0000000, 1.0000000, ], 100:[ 0.5000000, 0.5039794, 0.5397277, 0.6909132, 0.8401379, 0.9758939, 0.9999988, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, ], } probs = { 1: [ 0.500000000, 0.496817007, 0.468274483, 0.352416382, 0.250000000, 0.147583618, 0.062832958, 0.031725517, 0.015902251, 0.010606402, 0.006365349, 0.001591536, ], 10: [ 5.000000e-01, 4.961090e-01, 4.611604e-01, 3.139468e-01, 1.704466e-01, 3.669402e-02, 2.686668e-04, 7.947766e-07, 1.073031e-09, 1.980896e-11, 1.237155e-13, 1.200254e-19, ], 100:[ 5.000000e-01, 4.960206e-01, 4.602723e-01, 3.090868e-01, 1.598621e-01, 2.410609e-02, 1.225087e-06, 4.950844e-17, 4.997134e-37, 4.190166e-52, 7.236082e-73, 2.774197e-132, ], } for df in self.df: for x, p in zip(self.values, probs[df]): self.assertFloatEqualRel(t_high(x, df), p, eps=1e-4) for x, p in zip(self.negvalues, negprobs[df]): self.assertFloatEqualRel(t_high(x, df), p, eps=1e-4) def test_tprob(self): """tprob should match twice the t_high probability for abs(t)""" probs = { 1: [ 2*i for i in [ 0.500000000, 0.496817007, 0.468274483, 0.352416382, 0.250000000, 0.147583618, 0.062832958, 0.031725517, 0.015902251, 0.010606402, 0.006365349, 0.001591536, ]], 10: [ 2*i for i in [ 5.000000e-01, 4.961090e-01, 4.611604e-01, 3.139468e-01, 1.704466e-01, 3.669402e-02, 2.686668e-04, 7.947766e-07, 1.073031e-09, 1.980896e-11, 1.237155e-13, 1.200254e-19, ]], 100:[ 2*i for i in [ 5.000000e-01, 4.960206e-01, 4.602723e-01, 3.090868e-01, 1.598621e-01, 2.410609e-02, 1.225087e-06, 4.950844e-17, 4.997134e-37, 4.190166e-52, 7.236082e-73, 2.774197e-132, ]], } for df in self.df: for x, p in zip(self.values, probs[df]): self.assertFloatEqualRel(tprob(x, df), p, eps=1e-4) def test_poisson_low(self): """Lower tail of poisson should match R for integer successes""" #WARNING: Results only guaranteed for integer successes: floating #point _should_ yield reasonable values, but R rounds to int. expected = { (0, 0): 1, (0, 0.75): 0.4723666, (0, 1): 0.3678794, (0, 5): 0.006737947, (0, 113.7): 4.175586e-50, (2, 0): 1, (2, 3): 0.4231901, (2, 17.8): 3.296636e-06, (17, 29.6): 0.008753318, (180, 0): 1, (180, 137.4):0.999784, (180, 318):2.436995e-17, (180, 1024):8.266457e-233, } for (key, value) in expected.items(): self.assertFloatEqual(poisson_low(*key), value) def test_poisson_high(self): """Upper tail of poisson should match R for integer successes""" #WARNING: Results only guaranteed for integer successes: floating #point _should_ yield reasonable values, but R rounds to int. expected = { (0, 0): 0, (0, 0.75): 0.5276334, (0, 1): 0.6321206, (0, 5): 0.993262, (0, 113.7): 1, (2, 0): 0, (2, 3): 0.5768099, (2, 17.8): 0.9999967, (17, 29.6): 0.9912467, (180, 0): 0, (180, 137.4):0.0002159856, (180, 318):1, (180, 1024):1, } for (key, value) in expected.items(): self.assertFloatEqual(poisson_high(*key), value) def test_poisson_exact(self): """Poisson exact should match expected values from R""" expected = { (0, 0): 1, (0, 0.75): 0.4723666, (0, 1): 0.3678794, (0, 5): 0.006737947, (0, 113.7): 4.175586e-50, (2, 0): 0, (2, 3): 0.2240418, (2, 17.8): 2.946919e-06, (17, 29.6): 0.004034353, (180, 0): 0, (180, 137.4):7.287501e-05, (180, 318):1.067247e-17, (180, 1024):6.815085e-233, } for (key, value) in expected.items(): self.assertFloatEqual(poisson_exact(*key), value) def test_binomial_high(self): """Binomial high should match values from R for integer successes""" expected = { (0, 1, 0.5): 0.5, (1, 1, 0.5): 0, (1, 1, 0.0000001): 0, (1, 1, 0.9999999): 0, (3, 5, 0.75):0.6328125, (0, 60, 0.5): 1, (129, 130, 0.5):7.34684e-40, (299, 300, 0.099): 4.904089e-302, (9, 27, 0.0003): 4.958496e-29, (1032, 2050, 0.5): 0.3702155, (-1, 3, 0.1): 1, #if successes less than 0, return 1 (-0.5, 3, 0.1):1, } for (key, value) in expected.items(): self.assertFloatEqualRel(binomial_high(*key), value, 1e-4) #should reject if successes > trials or successes < -1 self.assertRaises(ValueError, binomial_high, 7, 5, 0.5) def test_binomial_low(self): """Binomial low should match values from R for integer successes""" expected = { (0, 1, 0.5): 0.5, (1, 1, 0.5): 1, (1, 1, 0.0000001): 1, (1, 1, 0.9999999): 1, (26, 50, .5): 0.6641, (3, 5, 0.75):0.3671875, (0, 60, 0.5): 8.673617e-19, (129, 130, 0.5):1, (299, 300, 0.099): 1, (9, 27, 0.0003): 1, (1032, 2050, 0.5): 0.6297845, } for (key, value) in expected.items(): self.assertFloatEqualRel(binomial_low(*key), value, 1e-4) def test_binomial_series(self): """binomial_exact should match values from R on a whole series""" expected = map(float, "0.0282475249 0.1210608210 0.2334744405 0.2668279320 0.2001209490 0.1029193452 0.0367569090 0.0090016920 0.0014467005 0.0001377810 0.0000059049".split()) for i in range(len(expected)): self.assertFloatEqual(binomial_exact(i, 10, 0.3), expected[i]) def test_binomial_exact(self): """binomial_exact should match values from R for integer successes""" expected = { (0, 1, 0.5): 0.5, (1, 1, 0.5): 0.5, (1, 1, 0.0000001): 1e-07, (1, 1, 0.9999999): 0.9999999, (3, 5, 0.75):0.2636719, (0, 60, 0.5): 8.673617e-19, (129, 130, 0.5):9.550892e-38, (299, 300, 0.099): 1.338965e-298, (9, 27, 0.0003): 9.175389e-26, (1032, 2050, 0.5): 0.01679804, } for (key, value) in expected.items(): self.assertFloatEqualRel(binomial_exact(*key), value, 1e-4) def test_binomial_exact_floats(self): """binomial_exact should be within limits for floating point numbers """ expected = { (18.3, 100, 0.2): (0.09089812, 0.09807429), (2.7,1050,0.006): (0.03615498, 0.07623827), (2.7,1050,0.06): (1.365299e-25, 3.044327e-24), (2,100.5,0.6): (7.303533e-37, 1.789727e-36), (10,100.5,.5):(7.578011e-18,1.365543e-17), (0.2, 60, 0.5): (8.673617e-19, 5.20417e-17), (.5,5,.3):(0.16807,0.36015), } for (key, value) in expected.items(): min_val, max_val = value assert min_val < binomial_exact(*key) < max_val #self.assertFloatEqualRel(binomial_exact(*key), value, 1e-4) def test_binomial_exact_errors(self): """binomial_exact should raise errors on invalid input""" self.assertRaises(ValueError, binomial_exact,10.2, 5, 0.33) self.assertRaises(ValueError, binomial_exact,-2, 5, 0.33) self.assertRaises(ValueError, binomial_exact, 10, 50, -2) self.assertRaises(ValueError, binomial_exact, 10, 50, 3) def test_f_high(self): """F high should match values from R for integer successes""" expected = { (1, 1, 0): 1, (1, 1, 1): 0.5, (1, 1, 20): 0.1400487, (1, 1, 1000000): 0.0006366196, (1, 10, 0): 1, (1,10, 5): 0.0493322, (1, 10, 20): 0.001193467, (10, 1, 0):1, (10, 10, 14.7): 0.0001062585, (13.7, 11.9, 3.8): 0.01340347, #test non-integer degrees of freedom #used following series to track down a bug after a failed test case (28, 29, 2): 0.03424088, (28, 29, 10): 1.053019e-08, (28, 29, 20): 1.628245e-12, (28, 29, 300): 5.038791e-29, (28, 35, 1): 0.4946777, (28, 37, 1): 0.4934486, (28, 38, 1): 0.4928721, (28, 38.001, 1): 0.4928716, (28, 38.5, 1): 0.4925927, (28, 39, 1): 0.492319, (28, 39, 10): 1.431901e-10, (28, 39, 20): 1.432014e-15, (28, 39, 30): 1.059964e-18, (28, 39, 50): 8.846678e-23, (28, 39, 10): 1.431901e-10, (28, 39, 300): 1.226935e-37, (28, 39, 50): 8.846678e-23, (28,39,304.7): 9.08154e-38, (28.4, 39.2, 304.7): 5.573927e-38, (1032, 2050, 0): 1, (1032, 2050, 4.15): 1.23535e-165, (1032, 2050, 0.5): 1, (1032, 2050, 0.1): 1, } e = expected.items() e.sort() for (key, value) in e: self.assertFloatEqualRel(f_high(*key), value) def test_f_low(self): """F low should match values from R for integer successes""" expected = { (1, 1, 0): 0, (1, 1, 1): 0.5, (1, 1, 20): 0.8599513, (1, 1, 1000000): 0.9993634, (1, 10, 0): 0, (1,10, 5): 0.9506678, (1, 10, 20): 0.9988065, (10, 1, 0):0, (10, 10, 14.7): 0.9998937, (28.4, 39.2, 304.7): 1, (1032, 2050, 0): 0, (1032, 2050, 4.15): 1, (1032, 2050, 0.5): 7.032663e-35, (1032, 2050, 0.1): 1.70204e-278, } for (key, value) in expected.items(): self.assertFloatEqualRel(f_low(*key), value) def test_fprob(self): """fprob should return twice the tail on a particular side""" error = 1e-4 #right-hand side self.assertFloatEqualAbs(fprob(10,10,1.2), 0.7788, eps=error) #left-hand side self.assertFloatEqualAbs(fprob(10,10,1.2, side='left'), 1.2212, eps=error) self.assertRaises(ValueError, fprob, 10,10,-3) self.assertRaises(ValueError, fprob, 10, 10, 1, 'non_valid_side') def test_stdtr(self): """stdtr should match cephes results""" t = [-10, -3.1, -0.5, -0.01, 0, 1, 0.5, 10] k = [2, 10, 100] exp = [ 0.00492622851166, 7.94776587798e-07, 4.9508444923e-17, 0.0451003650651, 0.00562532860804, 0.00125696358826, 0.333333333333, 0.313946802871, 0.309086782915, 0.496464554479, 0.496108987495, 0.496020605117, 0.5, 0.5, 0.5, 0.788675134595, 0.829553433849, 0.840137922108, 0.666666666667, 0.686053197129, 0.690913217085, 0.995073771488, 0.999999205223, 1.0, ] index = 0 for i in t: for j in k: self.assertFloatEqual(stdtr(j,i), exp[index]) index += 1 def test_bdtr(self): """bdtr should match cephes results""" k_s = [0,1,2,3,5] n_s = [5,10,1000] p_s = [1e-10, .1, .5, .9, .999999] exp = [ 0.9999999995, 0.59049, 0.03125, 1e-05, 1.00000000014e-30, 0.999999999, 0.3486784401, 0.0009765625, 1e-10, 1.00000000029e-60, 0.9999999, 1.74787125172e-46, 9.33263618503e-302, 0.0, 0.0, 1.0, 0.91854, 0.1875, 0.00046, 4.99999600058e-24, 1.0, 0.7360989291, 0.0107421875, 9.1e-09, 9.99999100259e-54, 1.0, 1.9595578811e-44, 9.34196882121e-299, 0.0, 0.0, 1.0, 0.99144, 0.5, 0.00856, 9.99998500087e-18, 1.0, 0.9298091736, 0.0546875, 3.736e-07, 4.49999200104e-47, 1.0, 1.09744951737e-42, 4.67099374325e-296, 0.0, 0.0, 1.0, 0.99954, 0.8125, 0.08146, 9.99998000059e-12, 1.0, 0.9872048016, 0.171875, 9.1216e-06, 1.19999685024e-40, 1.0, 4.09381247279e-41, 1.5554471507e-293, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9998530974, 0.623046875, 0.0016349374, 2.51998950038e-28, 1.0, 2.55654569306e-38, 7.7385053063e-289, 0.0, 0.0, ] index = 0 for k in k_s: for n in n_s: for p in p_s: self.assertFloatEqual(bdtr(k,n,p), exp[index]) index += 1 def test_bdtrc(self): """bdtrc should give same results as cephes""" k_s = [0,1,2,3,5] n_s = [5,10,1000] p_s = [1e-10, .1, .5, .9, .999999] exp = [ 4.999999999e-10, 0.40951, 0.96875, 0.99999, 1.0, 9.9999999955e-10, 0.6513215599, 0.9990234375, 0.9999999999, 1.0, 9.9999995005e-08, 1.0, 1.0, 1.0, 1.0, 9.999999998e-20, 0.08146, 0.8125, 0.99954, 1.0, 4.4999999976e-19, 0.2639010709, 0.9892578125, 0.9999999909, 1.0, 4.99499966766e-15, 1.0, 1.0, 1.0, 1.0, 9.9999999985e-30, 0.00856, 0.5, 0.99144, 1.0, 1.19999999937e-28, 0.0701908264, 0.9453125, 0.9999996264, 1.0, 1.66166987575e-22, 1.0, 1.0, 1.0, 1.0, 4.9999999996e-40, 0.00046, 0.1875, 0.91854, 0.99999999999, 2.09999999899e-38, 0.0127951984, 0.828125, 0.9999908784, 1.0, 4.14171214499e-30, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.09999999928e-58, 0.0001469026, 0.376953125, 0.9983650626, 1.0, 1.36817318242e-45, 1.0, 1.0, 1.0, 1.0, ] index = 0 for k in k_s: for n in n_s: for p in p_s: self.assertFloatEqual(bdtrc(k,n,p), exp[index]) index += 1 def test_pdtr(self): """pdtr should match cephes results""" k_s = [0,1,2,5,10] m_s = [1e-9, 0.1,0.5,1,2,31] exp = [ 0.999999999 , 0.904837418036 , 0.606530659713 , 0.367879441171 , 0.135335283237 , 3.44247710847e-14 , 1.0 , 0.99532115984 , 0.909795989569 , 0.735758882343 , 0.40600584971 , 1.10159267471e-12 , 1.0 , 0.99984534693 , 0.985612322033 , 0.919698602929 , 0.676676416183 , 1.76426951809e-11 , 1.0 , 0.999999998725 , 0.999985835063 , 0.999405815182 , 0.983436391519 , 9.72616712615e-09 , 1.0 , 1.0 , 0.999999999992 , 0.999999989952 , 0.999991691776 , 1.12519146046e-05 , ] index = 0 for k in k_s: for m in m_s: self.assertFloatEqual(pdtr(k,m), exp[index]) index += 1 def test_pdtrc(self): """pdtrc should match cephes results""" k_s = [0,1,2,5,10] m_s = [1e-9, 0.1,0.5,1,2,31] exp = [ 9.999999995e-10 , 0.095162581964 , 0.393469340287 , 0.632120558829 , 0.864664716763 , 1.0 , 4.99999999667e-19 , 0.00467884016044 , 0.090204010431 , 0.264241117657 , 0.59399415029 , 0.999999999999 , 1.66666666542e-28 , 0.000154653070265 , 0.014387677967 , 0.0803013970714 , 0.323323583817 , 0.999999999982 , 1.3888888877e-57 , 1.27489869223e-09 , 1.41649373223e-05 , 0.000594184817582 , 0.0165636084806 , 0.999999990274 , 2.50521083625e-107 , 2.28584493079e-19 , 7.74084073923e-12 , 1.00477663757e-08 , 8.30822436848e-06 , 0.999988748085 , ] index = 0 for k in k_s: for m in m_s: self.assertFloatEqual(pdtrc(k,m), exp[index]) index += 1 def test_fdtr(self): """fdtr should match cephes results""" a_s = [1, 2, 10, 1000] b_s = a_s x_s = [0, 0.01, 0.5, 10, 521.4] exp = [ 0.0, 0.0634510348611, 0.391826552031, 0.805017770958, 0.972137685271, 0.0, 0.0705345615859, 0.4472135955, 0.912870929175, 0.998087586699, 0.0, 0.0776792814356, 0.504352495617, 0.989880440265, 0.999999999415, 0.0, 0.0796356309764, 0.520335137562, 0.998387447605, 1.0, 0.0, 0.00985245702333, 0.292893218813, 0.781782109764, 0.96904781206, 0.0, 0.00990099009901, 0.333333333333, 0.909090909091, 0.99808575804, 0.0, 0.00994027888402, 0.379078676941, 0.995884773663, 0.999999999923, 0.0, 0.00995006724716, 0.393317789705, 0.999949891187, 1.0, 0.0, 1.5895531756e-06, 0.18766987087, 0.758331535711, 0.965930763936, 0.0, 2.44851927021e-07, 0.185934432082, 0.90573080983, 0.998084291751, 0.0, 1.15978163168e-08, 0.144845806026, 0.999428447457, 0.999999999997, 0.0, 2.54720538101e-09, 0.109321108726, 1.0, 1.0, 0.0, 1.66707029586e-22, 0.157610464133, 0.751895627261, 0.965077362955, 0.0, 2.56671102571e-40, 0.135876263477, 0.904846465249, 0.998083928382, 0.0, 5.1131107462e-143, 0.030392376141, 0.999825108037, 0.999999999999, 0.0, 0.0, 9.96002853277e-28, 1.0, 1.0, ] index = 0 for a in a_s: for b in b_s: for x in x_s: self.assertFloatEqual(fdtr(a,b,x), exp[index]) index += 1 def test_fdtrc(self): """fdtrc should match cephes results""" a_s = [1, 2, 10, 1000] b_s = a_s x_s = [0, 0.01, 0.5, 10, 521.4] exp = [ 1.0, 0.936548965139, 0.608173447969, 0.194982229042, 0.0278623147287, 1.0, 0.929465438414, 0.5527864045, 0.0871290708247, 0.00191241330122, 1.0, 0.922320718564, 0.495647504383, 0.0101195597354, 5.85364343244e-10, 1.0, 0.920364369024, 0.479664862438, 0.00161255239482, 3.24963344513e-93, 1.0, 0.990147542977, 0.707106781187, 0.218217890236, 0.0309521879405, 1.0, 0.990099009901, 0.666666666667, 0.0909090909091, 0.00191424196018, 1.0, 0.990059721116, 0.620921323059, 0.00411522633745, 7.73162209771e-11, 1.0, 0.990049932753, 0.606682210295, 5.01088134545e-05, 7.71037335669e-156, 1.0, 0.999998410447, 0.81233012913, 0.241668464289, 0.0340692360638, 1.0, 0.999999755148, 0.814065567918, 0.0942691901701, 0.00191570824928, 1.0, 0.999999988402, 0.855154193974, 0.000571552543402, 3.21796660031e-12, 1.0, 0.999999997453, 0.890678891274, 3.96065609687e-16, 0.0, 1.0, 1.0, 0.842389535866, 0.248104372739, 0.0349226370457, 1.0, 1.0, 0.864123736523, 0.0951535347509, 0.00191607161849, 1.0, 1.0, 0.969607623859, 0.00017489196271, 6.83862415869e-13, 1.0, 1.0, 1.0, 6.68418402018e-243, 0.0, ] index = 0 for a in a_s: for b in b_s: for x in x_s: self.assertFloatEqual(fdtrc(a,b,x), exp[index]) index += 1 def test_gdtr(self): """gdtr should match cephes results""" a_s = [1, 2, 10, 1000] b_s = a_s x_s = [0, 0.01, 0.5, 10, 521.4] exp = [ 0.0, 0.00995016625083, 0.393469340287, 0.99995460007, 1.0, 0.0, 4.96679133403e-05, 0.090204010431, 0.999500600773, 1.0, 0.0, 2.7307942837e-27, 1.70967002935e-10, 0.542070285528, 1.0, 0.0, 0.0, 0.0, 0.0, 2.78154480191e-77, 0.0, 0.0198013266932, 0.632120558829, 0.999999997939, 1.0, 0.0, 0.00019735322711, 0.264241117657, 0.999999956716, 1.0, 0.0, 2.77103020131e-24, 1.11425478339e-07, 0.995004587692, 1.0, 0.0, 0.0, 0.0, 0.0, 0.91070640569, 0.0, 0.095162581964, 0.993262053001, 1.0, 1.0, 0.0, 0.00467884016044, 0.959572318005, 1.0, 1.0, 0.0, 2.51634780677e-17, 0.0318280573062, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.99995460007, 1.0, 1.0, 1.0, 0.0, 0.999500600773, 1.0, 1.0, 1.0, 0.0, 0.542070285528, 1.0, 1.0, 1.0, 0.0, 0.0, 3.29827279707e-86, 1.0, 1.0, ] index = 0 for a in a_s: for b in b_s: for x in x_s: self.assertFloatEqual(gdtr(a,b,x), exp[index]) index += 1 def test_gdtrc(self): """gdtrc should match cephes results""" a_s = [1, 2, 10, 1000] b_s = a_s x_s = [0, 0.01, 0.5, 10, 521.4] exp = [ 1.0, 0.990049833749, 0.606530659713, 4.53999297625e-05, 3.62123855523e-227, 1.0, 0.999950332087, 0.909795989569, 0.000499399227387, 1.89173502125e-224, 1.0, 1.0, 0.999999999829, 0.457929714472, 2.89188102723e-208, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.980198673307, 0.367879441171, 2.06115362244e-09, 0.0, 1.0, 0.999802646773, 0.735758882343, 4.32842260712e-08, 0.0, 1.0, 1.0, 0.999999888575, 0.00499541230831, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0892935943104, 1.0, 0.904837418036, 0.00673794699909, 3.72007597602e-44, 0.0, 1.0, 0.99532115984, 0.0404276819945, 3.75727673578e-42, 0.0, 1.0, 1.0, 0.968171942694, 1.12534739608e-31, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 4.53999297625e-05, 7.12457640674e-218, 0.0, 0.0, 1.0, 0.000499399227387, 3.56941277978e-215, 0.0, 0.0, 1.0, 0.457929714472, 3.90479663912e-199, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, ] index = 0 for a in a_s: for b in b_s: for x in x_s: self.assertFloatEqual(gdtrc(a,b,x), exp[index]) index += 1 def test_chdtri(self): """chdtri should match cephes results""" k_s = [1,2,5,10,100] p_s = [1e-50, 1e-9, .02, .5, .8, .99] exp = [ 224.384748319, 37.3248930514, 5.41189443105, 0.45493642312, 0.0641847546673, 0.00015708785791, 230.258509299, 41.4465316739, 7.82404601086, 1.38629436112, 0.446287102628, 0.020100671707, 244.127298027, 50.6921937015, 13.388222599, 4.3514601911, 2.34253430584, 0.554298076728, 262.995620961, 62.9454574206, 21.1607675413, 9.34181776559, 6.17907925604, 2.55821216019, 478.347499744, 209.317598707, 131.141676866, 99.334129236, 87.9453359228, 70.0648949254, ] index = 0 for k in k_s: for p in p_s: self.assertFloatEqual(chdtri(k,p), exp[index]) index += 1 def test_stdtri(self): """stdtri should match cephes results""" k_s = [1,2,5,10,100] p_s = [1e-50, 1e-9, .02, .5, .8, .99] exp = [ -3.18309886184e+49, -318309886.184, -15.8945448441, 8.1775627727e-17, 1.37638192049, 31.8205159538, -7.07106781216e+24, -22360.6797414, -4.84873221444, 7.48293180888e-17, 1.06066017178, 6.96455671876, -15683925591.1, -98.9372246484, -2.75650852191, 6.976003623e-17, 0.919543780236, 3.36492999891, -256452.571877, -20.1446977667, -2.35931462368, 6.80574793291e-17, 0.879057828551, 2.76376945745, -28.9584072963, -6.59893982023, -2.08088390123, 6.6546053747e-17, 0.845230424487, 2.3642173659, ] index = 0 for k in k_s: for p in p_s: self.assertFloatEqual(stdtri(k,p), exp[index]) index += 1 def test_pdtri(self): """pdtri should match cephes results""" k_s = [1,2,5,10,100] p_s = [1e-50, 1e-9, .02, .5, .8, .99] exp = [ 119.924420375, 23.9397278656, 5.83392170192, 1.67834699002, 0.824388309033, 0.148554740253, 124.094307191, 26.6722865587, 7.51660387561, 2.67406031372, 1.53504420264, 0.436045165078, 134.901981814, 33.6746016741, 12.0269783451, 5.67016118871, 3.90366383933, 1.7852844853, 150.2138305, 43.627975401, 18.8297496417, 10.6685224038, 8.15701989758, 4.77124616939, 332.371212972, 173.368244558, 122.695978128, 100.666862949, 92.4593447729, 79.0999186597, ] index = 0 for k in k_s: for p in p_s: self.assertFloatEqual(pdtri(k,p), exp[index]) index += 1 def test_bdtri(self): """bdtri should match cephes results""" k_s = [0,1,2,3] n_s = [5,10,1000] p_s = [1e-10, .1, .5, .9, .999999] exp = [ 0.99, 0.36904265552, 0.129449436704, 0.020851637639, 2.00000080006e-07, 0.9, 0.205671765276, 0.0669670084632, 0.0104807417938, 1.00000045003e-07, 0.0227627790442, 0.00229993617745, 0.000692907009547, 0.000105354965434, 1.00000049953e-09, 0.997884361719, 0.58389037462, 0.313810170456, 0.112234958546, 0.000316327821398, 0.939678616058, 0.336847723307, 0.162262728195, 0.0545286199977, 0.00014913049349, 0.0260030189545, 0.0038841043984, 0.00167777786542, 0.000531936197341, 1.415587631e-06, 0.999784533318, 0.753363546712, 0.5, 0.246636453288, 0.00465241636163, 0.964779035441, 0.449603888674, 0.258574723285, 0.11582527803, 0.00203463563411, 0.0287538329681, 0.00531348536403, 0.00267315927217, 0.00110256069953, 1.82723947076e-05, 0.999996837712, 0.887765041454, 0.686189829544, 0.41610962538, 0.0212382182007, 0.981054188003, 0.551730832384, 0.355099967912, 0.187562296647, 0.00839131408953, 0.0312483560212, 0.00666849533707, 0.00367082709364, 0.00174586632568, 7.10965576424e-05, ] index = 0 for k in k_s: for n in n_s: for p in p_s: self.assertFloatEqual(bdtri(k,n,p), exp[index]) index += 1 def test_gdtri(self): """gdtri should match cephes results""" k_s = [1,2,4,10,100] n_s = k_s p_s = [1e-9, .02, .5, .8, .99] exp = [ 1.0000000005e-09, 0.0202027073175, 0.69314718056, 1.60943791243, 4.60517018599, 4.47220262303e-05, 0.214699095008, 1.67834699002, 2.994308347, 6.63835206799, 0.0124777531242, 1.01623845904, 3.67206074885, 5.51504571515, 10.0451175148, 0.602134838869, 4.61834927756, 9.66871461471, 12.5187528198, 18.7831173933, 51.1433022288, 80.5501391278, 99.6668649193, 108.304391619, 124.722561491, 5.0000000025e-10, 0.0101013536588, 0.34657359028, 0.804718956217, 2.30258509299, 2.23610131152e-05, 0.107349547504, 0.839173495008, 1.4971541735, 3.319176034, 0.00623887656209, 0.50811922952, 1.83603037443, 2.75752285758, 5.02255875742, 0.301067419435, 2.30917463878, 4.83435730736, 6.25937640991, 9.39155869666, 25.5716511144, 40.2750695639, 49.8334324597, 54.1521958095, 62.3612807454, 2.50000000125e-10, 0.00505067682938, 0.17328679514, 0.402359478109, 1.1512925465, 1.11805065576e-05, 0.053674773752, 0.419586747504, 0.748577086751, 1.659588017, 0.00311943828105, 0.25405961476, 0.918015187213, 1.37876142879, 2.51127937871, 0.150533709717, 1.15458731939, 2.41717865368, 3.12968820495, 4.69577934833, 12.7858255572, 20.1375347819, 24.9167162298, 27.0760979048, 31.1806403727, 1.0000000005e-10, 0.00202027073175, 0.069314718056, 0.160943791243, 0.460517018599, 4.47220262303e-06, 0.0214699095008, 0.167834699002, 0.2994308347, 0.663835206799, 0.00124777531242, 0.101623845904, 0.367206074885, 0.551504571515, 1.00451175148, 0.0602134838869, 0.461834927756, 0.966871461471, 1.25187528198, 1.87831173933, 5.11433022288, 8.05501391278, 9.96668649193, 10.8304391619, 12.4722561491, 1.0000000005e-11, 0.000202027073175, 0.0069314718056, 0.0160943791243, 0.0460517018599, 4.47220262303e-07, 0.00214699095008, 0.0167834699002, 0.02994308347, 0.0663835206799, 0.000124777531242, 0.0101623845904, 0.0367206074885, 0.0551504571515, 0.100451175148, 0.00602134838869, 0.0461834927756, 0.0966871461471, 0.125187528198, 0.187831173933, 0.511433022288, 0.805501391278, 0.996668649193, 1.08304391619, 1.24722561491, ] index = 0 for k in k_s: for n in n_s: for p in p_s: self.assertFloatEqual(gdtri(k,n,p), exp[index]) index += 1 def test_fdtri(self): """fdtri should match cephes results""" k_s = [1,2,4,10,100] n_s = k_s p_s = [1e-50, 1e-9, .02, .5, .8, .99] exp = [ 0.0, 2.46740096071e-18, 0.000987610197427, 1.0, 9.472135955, 4052.18069548, 0.0, 1.99999988687e-18, 0.000800320128051, 0.666666666667, 3.55555555556, 98.5025125628, 0.0, 1.77777767722e-18, 0.000711321880645, 0.548632170413, 2.35072147881, 21.1976895844, 0.0, 1.65119668161e-18, 0.000660638708985, 0.489736921158, 1.88288794493, 10.0442892734, 0.0, 1.57866975531e-18, 0.000631602221127, 0.458262714634, 1.66429288986, 6.89530103058, 0.0, 9.99999973218e-10, 0.0206164098292, 1.5, 12.0, 4999.5, 0.0, 9.99999972718e-10, 0.0204081632653, 1.0, 4.0, 99.0, 0.0, 9.99999972468e-10, 0.0203050891044, 0.828427124746, 2.472135955, 18.0, 0.0, 9.99999972318e-10, 0.0202435772829, 0.743491774985, 1.89864830731, 7.55943215755, 0.0, 9.99999972228e-10, 0.0202067893611, 0.697973989501, 1.63562099482, 4.82390980716, 0.0, 1.29104998825e-05, 0.0712270257663, 1.82271484235, 13.6443218387, 5624.58332963, 0.0, 1.58118880931e-05, 0.082357834815, 1.20710678119, 4.2360679775, 99.2493718553, 0.0, 1.82578627816e-05, 0.0917479893415, 1.0, 2.48261291932, 15.9770248526, 0.0, 2.04128031324e-05, 0.0999726146531, 0.898817134423, 1.82861100515, 5.99433866163, 0.0, 2.21407117017e-05, 0.106518545067, 0.844891468084, 1.5273126184, 3.5126840636, 0.0, 0.00213897888638, 0.130917099116, 2.04191262042, 14.7718897826, 6055.8467074, 0.0, 0.00322083313175, 0.168531162323, 1.34500479177, 4.38216390487, 99.3991959745, 0.0, 0.00448830777955, 0.207656634378, 1.11257336081, 2.45957986729, 14.5459008033, 0.0, 0.00608578074458, 0.251574092492, 1.0, 1.73159473193, 4.84914680208, 0.0, 0.00800159033308, 0.298648905106, 0.940477156977, 1.38089597558, 2.50331112688, 0.0, 3.09672866088e-11, 0.178906118636, 2.18215440197, 15.4973240414, 6334.110036, 0.0, 5.2776234633e-11, 0.2457526061, 1.43271814572, 4.47142755584, 99.4891628084, 0.0, 0.0659164713677, 0.326585865322, 1.18358397235, 2.43020291912, 13.5769915067, 0.0, 0.119865858243, 0.442669184276, 1.06329004653, 1.63265061785, 4.01371941549, 0.0, 0.289673110482, 0.661509869668, 1.0, 1.1839371445, 1.59766912303, ] index = 0 for k in k_s: for n in n_s: for p in p_s: self.assertFloatEqual(fdtri(k,n,p), exp[index]) index += 1 if __name__ == "__main__": main()