#!/usr/bin/env python """Unit tests for statistical tests and utility functions. """ from cogent.util.unit_test import TestCase, main from cogent.maths.stats.test import tail, G_2_by_2,G_fit, likelihoods,\ posteriors, bayes_updates, t_paired, t_one_sample, t_two_sample, \ t_one_observation,correlation, correlation_matrix, z_test, z_tailed_prob, \ t_tailed_prob, sign_test,\ reverse_tails, ZeroExpectedError, combinations, multiple_comparisons, \ multiple_inverse, multiple_n, fisher, regress, regress_major,\ f_value, f_two_sample, calc_contingency_expected, G_fit_from_Dict2D, \ chi_square_from_Dict2D, MonteCarloP, \ regress_residuals, safe_sum_p_log_p, G_ind, regress_origin, stdev_from_mean, \ regress_R2, permute_2d, mantel, kendall_correlation, std, median,\ get_values_from_matrix, get_ltm_cells, distance_matrix_permutation_test,\ ANOVA_one_way from numpy import array, reshape, arange, ones, testing, cov, sqrt from cogent.util.dict2d import Dict2D import math from cogent.maths.stats.util import Numbers __author__ = "Rob Knight" __copyright__ = "Copyright 2007-2009, The Cogent Project" __credits__ = ["Rob Knight", "Catherine Lozupone", "Gavin Huttley", "Sandra Smit", "Daniel McDonald"] __license__ = "GPL" __version__ = "1.4.1" __maintainer__ = "Rob Knight" __email__ = "rob@spot.colorado.edu" __status__ = "Production" class TestsTests(TestCase): """Tests miscellaneous functions.""" def test_std(self): """Should produce a standard deviation of 1.0 for a std normal dist""" expected = 1.58113883008 self.assertFloatEqual(std(array([1,2,3,4,5])), expected) expected_a = array([expected, expected, expected, expected, expected]) a = array([[1,2,3,4,5],[5,1,2,3,4],[4,5,1,2,3],[3,4,5,1,2],[2,3,4,5,1]]) self.assertFloatEqual(std(a,axis=0), expected_a) self.assertFloatEqual(std(a,axis=1), expected_a) self.assertRaises(ValueError, std, a, 5) def test_std_2d(self): """Should produce from 2darray the same stdevs as scipy.stats.std""" inp = array([[1,2,3],[4,5,6]]) exps = ( #tuple(scipy_std(inp, ax) for ax in [None, 0, 1]) 1.8708286933869707, array([ 2.12132034, 2.12132034, 2.12132034]), array([ 1., 1.])) results = tuple(std(inp, ax) for ax in [None, 0, 1]) for obs, exp in zip(results, exps): testing.assert_almost_equal(obs, exp) def test_std_3d(self): """Should produce from 3darray the same std devs as scipy.stats.std""" inp3d = array(#2,2,3 [[[ 0, 2, 2], [ 3, 4, 5]], [[ 1, 9, 0], [ 9, 10, 1]]]) exp3d = (#for axis None, 0, 1, 2: calc from scipy.stats.std 3.63901418552, array([[ 0.70710678, 4.94974747, 1.41421356], [ 4.24264069, 4.24264069, 2.82842712]]), array([[ 2.12132034, 1.41421356, 2.12132034], [ 5.65685425, 0.70710678, 0.70710678]]), array([[ 1.15470054, 1. ], [ 4.93288286, 4.93288286]])) res = tuple(std(inp3d, ax) for ax in [None, 0, 1, 2]) for obs, exp in zip(res, exp3d): testing.assert_almost_equal(obs, exp) def test_median(self): """_median should work similarly to numpy.mean (in terms of axis)""" m = array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]) expected = 6.5 observed = median(m, axis=None) self.assertEqual(observed, expected) expected = array([5.5, 6.5, 7.5]) observed = median(m, axis=0) self.assertEqual(observed, expected) expected = array([2.0, 5.0, 8.0, 11.0]) observed = median(m, axis=1) self.assertEqual(observed, expected) self.assertRaises(ValueError, median, m, 10) def test_tail(self): """tail should return x/2 if test is true; 1-(x/2) otherwise""" self.assertFloatEqual(tail(0.25, 'a'=='a'), 0.25/2) self.assertFloatEqual(tail(0.25, 'a'!='a'), 1-(0.25/2)) def test_combinations(self): """combinations should return correct binomial coefficient""" self.assertFloatEqual(combinations(5,3), 10) self.assertFloatEqual(combinations(5,2), 10) #only one way to pick no items or the same number of items self.assertFloatEqual(combinations(123456789, 0), 1) self.assertFloatEqual(combinations(123456789, 123456789), 1) #n ways to pick one item self.assertFloatEqual(combinations(123456789, 1), 123456789) #n(n-1)/2 ways to pick 2 items self.assertFloatEqual(combinations(123456789, 2), 123456789*123456788/2) #check an arbitrary value in R self.assertFloatEqual(combinations(1234567, 12), 2.617073e64) def test_multiple_comparisons(self): """multiple_comparisons should match values from R""" self.assertFloatEqual(multiple_comparisons(1e-7, 10000), 1-0.9990005) self.assertFloatEqual(multiple_comparisons(0.05, 10), 0.4012631) self.assertFloatEqual(multiple_comparisons(1e-20, 1), 1e-20) self.assertFloatEqual(multiple_comparisons(1e-300, 1), 1e-300) self.assertFloatEqual(multiple_comparisons(0.95, 3),0.99987499999999996) self.assertFloatEqual(multiple_comparisons(0.75, 100),0.999999999999679) self.assertFloatEqual(multiple_comparisons(0.5, 1000),1) self.assertFloatEqual(multiple_comparisons(0.01, 1000),0.99995682875259) self.assertFloatEqual(multiple_comparisons(0.5, 5), 0.96875) self.assertFloatEqual(multiple_comparisons(1e-20, 10), 1e-19) def test_multiple_inverse(self): """multiple_inverse should invert multiple_comparisons results""" #NOTE: multiple_inverse not very accurate close to 1 self.assertFloatEqual(multiple_inverse(1-0.9990005, 10000), 1e-7) self.assertFloatEqual(multiple_inverse(0.4012631 , 10), 0.05) self.assertFloatEqual(multiple_inverse(1e-20, 1), 1e-20) self.assertFloatEqual(multiple_inverse(1e-300, 1), 1e-300) self.assertFloatEqual(multiple_inverse(0.96875, 5), 0.5) self.assertFloatEqual(multiple_inverse(1e-19, 10), 1e-20) def test_multiple_n(self): """multiple_n should swap parameters in multiple_comparisons""" self.assertFloatEqual(multiple_n(1e-7, 1-0.9990005), 10000) self.assertFloatEqual(multiple_n(0.05, 0.4012631), 10) self.assertFloatEqual(multiple_n(1e-20, 1e-20), 1) self.assertFloatEqual(multiple_n(1e-300, 1e-300), 1) self.assertFloatEqual(multiple_n(0.95,0.99987499999999996),3) self.assertFloatEqual(multiple_n(0.5,0.96875),5) self.assertFloatEqual(multiple_n(1e-20, 1e-19), 10) def test_fisher(self): """fisher results should match p 795 Sokal and Rohlf""" self.assertFloatEqual(fisher([0.073,0.086,0.10,0.080,0.060]), 0.0045957946540917905) def test_regress(self): """regression slope, intercept should match p 459 Sokal and Rohlf""" x = [0, 12, 29.5,43,53,62.5,75.5,85,93] y = [8.98, 8.14, 6.67, 6.08, 5.90, 5.83, 4.68, 4.20, 3.72] self.assertFloatEqual(regress(x, y), (-0.05322, 8.7038), 0.001) #higher precision from OpenOffice self.assertFloatEqual(regress(x, y), (-0.05322215,8.70402730)) def test_regress_origin(self): """regression slope constrained through origin should match Excel""" x = array([1,2,3,4]) y = array([4,2,6,8]) self.assertFloatEqual(regress_origin(x, y), (1.9333333,0)) def test_regress_R2(self): """regress_R2 returns the R^2 value of a regression""" x = [1.0,2.0,3.0,4.0,5.0] y = [2.1,4.2,5.9,8.4,9.6] result = regress_R2(x, y) self.assertFloatEqual(result, 0.99171419347896) def test_regress_residuals(self): """regress_residuals reprts error for points in linear regression""" x = [1.0,2.0,3.0,4.0,5.0] y = [2.1,4.2,5.9,8.4,9.6] result = regress_residuals(x, y) self.assertFloatEqual(result, [-0.1, 0.08, -0.14, 0.44, -0.28]) def test_stdev_from_mean(self): """stdev_from_mean returns num std devs from mean for each val in x""" x = [2.1, 4.2, 5.9, 8.4, 9.6] result = stdev_from_mean(x) self.assertFloatEqual(result, [-1.292463399014413, -0.60358696806764478, -0.045925095396451399, 0.77416589382589174, 1.1678095686526162]) def test_regress_major(self): """major axis regression should match p 589 Sokal and Rohlf""" #Note that the Sokal and Rohlf example flips the axes, such that the #equation is for explaining x in terms of y, not y in terms of x. #Behavior here is the reverse, for easy comparison with regress. y = [159, 179, 100, 45, 384, 230, 100, 320, 80, 220, 320, 210] x = [14.40, 15.20, 11.30, 2.50, 22.70, 14.90, 1.41, 15.81, 4.19, 15.39, 17.25, 9.52] self.assertFloatEqual(regress_major(x, y), (18.93633,-32.55208)) def test_sign_test(self): """sign_test, should match values from R""" v = [("two sided", 26, 50, 0.88772482734078251), ("less", 26, 50, 0.6641), ("l", 10, 50, 1.193066583837777e-05), ("hi", 30, 50, 0.1013193755322703), ("h", 0, 50, 1.0), ("2", 30, 50, 0.20263875106454063), ("h", 49, 50, 4.5297099404706387e-14), ("h", 50, 50, 8.8817841970012543e-16) ] for alt, success, trials, p in v: result = sign_test(success, trials, alt=alt) self.assertFloatEqual(result, p, eps=1e-5) class GTests(TestCase): """Tests implementation of the G tests for fit and independence.""" def test_G_2_by_2_2tailed_equal(self): """G_2_by_2 should return 0 if all cell counts are equal""" self.assertFloatEqual(0, G_2_by_2(1, 1, 1, 1, False, False)[0]) self.assertFloatEqual(0, G_2_by_2(100, 100, 100, 100, False, False)[0]) self.assertFloatEqual(0, G_2_by_2(100, 100, 100, 100, True, False)[0]) def test_G_2_by_2_bad_data(self): """G_2_by_2 should raise ValueError if any counts are negative""" self.assertRaises(ValueError, G_2_by_2, 1, -1, 1, 1) def test_G_2_by_2_2tailed_examples(self): """G_2_by_2 values should match examples in Sokal & Rohlf""" #example from p 731, Sokal and Rohlf (1995) #without correction self.assertFloatEqual(G_2_by_2(12, 22, 16, 50, False, False)[0], 1.33249, 0.0001) self.assertFloatEqual(G_2_by_2(12, 22, 16, 50, False, False)[1], 0.24836, 0.0001) #with correction self.assertFloatEqual(G_2_by_2(12, 22, 16, 50, True, False)[0], 1.30277, 0.0001) self.assertFloatEqual(G_2_by_2(12, 22, 16, 50, True, False)[1], 0.25371, 0.0001) def test_G_2_by_2_1tailed_examples(self): """G_2_by_2 values should match values from codon_binding program""" #first up...the famous arginine case self.assertFloatEqualAbs(G_2_by_2(36, 16, 38, 106), (29.111609, 0), 0.00001) #then some other miscellaneous positive and negative values self.assertFloatEqualAbs(G_2_by_2(0,52,12,132), (-7.259930, 0.996474), 0.00001) self.assertFloatEqualAbs(G_2_by_2(5,47,14,130), (-0.000481, 0.508751), 0.00001) self.assertFloatEqualAbs(G_2_by_2(5,47,36,108), (-6.065167, 0.993106), 0.00001) def test_calc_contingency_expected(self): """calcContingencyExpected returns new matrix with expected freqs""" matrix = Dict2D({'rest_of_tree': {'env1': 2, 'env3': 1, 'env2': 0}, 'b': {'env1': 1, 'env3': 1, 'env2': 3}}) result = calc_contingency_expected(matrix) self.assertFloatEqual(result['rest_of_tree']['env1'], [2, 1.125]) self.assertFloatEqual(result['rest_of_tree']['env3'], [1, 0.75]) self.assertFloatEqual(result['rest_of_tree']['env2'], [0, 1.125]) self.assertFloatEqual(result['b']['env1'], [1, 1.875]) self.assertFloatEqual(result['b']['env3'], [1, 1.25]) self.assertFloatEqual(result['b']['env2'], [3, 1.875]) def test_Gfit_unequal_lists(self): """Gfit should raise errors if lists unequal""" #lists must be equal self.assertRaises(ValueError, G_fit, [1, 2, 3], [1, 2]) def test_Gfit_negative_observeds(self): """Gfit should raise ValueError if any observeds are negative.""" self.assertRaises(ValueError, G_fit, [-1, 2, 3], [1, 2, 3]) def test_Gfit_nonpositive_expecteds(self): """Gfit should raise ZeroExpectedError if expecteds are zero/negative""" self.assertRaises(ZeroExpectedError, G_fit, [1, 2, 3], [0, 1, 2]) self.assertRaises(ZeroExpectedError, G_fit, [1, 2, 3], [-1, 1, 2]) def test_Gfit_good_data(self): """Gfit tests for fit should match examples in Sokal and Rohlf""" #example from p. 699, Sokal and Rohlf (1995) obs = [63, 31, 28, 12, 39, 16, 40, 12] exp = [ 67.78125, 22.59375, 22.59375, 7.53125, 45.18750, 15.06250, 45.18750, 15.06250] #without correction self.assertFloatEqualAbs(G_fit(obs, exp, False)[0], 8.82397, 0.00002) self.assertFloatEqualAbs(G_fit(obs, exp, False)[1], 0.26554, 0.00002) #with correction self.assertFloatEqualAbs(G_fit(obs, exp)[0], 8.76938, 0.00002) self.assertFloatEqualAbs(G_fit(obs, exp)[1], 0.26964, 0.00002) #example from p. 700, Sokal and Rohlf (1995) obs = [130, 46] exp = [132, 44] #without correction self.assertFloatEqualAbs(G_fit(obs, exp, False)[0], 0.12002, 0.00002) self.assertFloatEqualAbs(G_fit(obs, exp, False)[1], 0.72901, 0.00002) #with correction self.assertFloatEqualAbs(G_fit(obs, exp)[0], 0.11968, 0.00002) self.assertFloatEqualAbs(G_fit(obs, exp)[1], 0.72938, 0.00002) def test_safe_sum_p_log_p(self): """safe_sum_p_log_p should ignore zero elements, not raise error""" m = array([2,4,0,8]) self.assertEqual(safe_sum_p_log_p(m,2), 2*1+4*2+8*3) def test_G_ind(self): """G test for independence should match Sokal and Rohlf p 738 values""" a = array([[29,11],[273,191],[8,31],[64,64]]) self.assertFloatEqual(G_ind(a)[0], 28.59642) self.assertFloatEqual(G_ind(a, True)[0], 28.31244) def test_G_fit_from_Dict2D(self): """G_fit_from_Dict2D runs G-fit on data in a Dict2D """ matrix = Dict2D({'Marl': {'val':[2, 5.2]}, 'Chalk': {'val':[10, 5.2]}, 'Sandstone':{'val':[8, 5.2]}, 'Clay':{'val':[2, 5.2]}, 'Limestone':{'val':[4, 5.2]} }) g_val, prob = G_fit_from_Dict2D(matrix) self.assertFloatEqual(g_val, 9.84923) self.assertFloatEqual(prob, 0.04304536) def test_chi_square_from_Dict2D(self): """chi_square_from_Dict2D calcs a Chi-Square and p value from Dict2D""" #test1 obs_matrix = Dict2D({'rest_of_tree': {'env1': 2, 'env3': 1, 'env2': 0}, 'b': {'env1': 1, 'env3': 1, 'env2': 3}}) input_matrix = calc_contingency_expected(obs_matrix) test, csp = chi_square_from_Dict2D(input_matrix) self.assertFloatEqual(test, 3.0222222222222221) #test2 test_matrix_2 = Dict2D({'Marl': {'val':[2, 5.2]}, 'Chalk': {'val':[10, 5.2]}, 'Sandstone':{'val':[8, 5.2]}, 'Clay':{'val':[2, 5.2]}, 'Limestone':{'val':[4, 5.2]} }) test2, csp2 = chi_square_from_Dict2D(test_matrix_2) self.assertFloatEqual(test2, 10.1538461538) self.assertFloatEqual(csp2, 0.0379143890013) #test3 matrix3_obs = Dict2D({'AIDS':{'Males':4, 'Females':2, 'Both':3}, 'No_AIDS':{'Males':3, 'Females':16, 'Both':2} }) matrix3 = calc_contingency_expected(matrix3_obs) test3, csp3 = chi_square_from_Dict2D(matrix3) self.assertFloatEqual(test3, 7.6568405139833722) self.assertFloatEqual(csp3, 0.0217439383468) class LikelihoodTests(TestCase): """Tests implementations of likelihood calculations.""" def test_likelihoods_unequal_list_lengths(self): """likelihoods should raise ValueError if input lists unequal length""" self.assertRaises(ValueError, likelihoods, [1, 2], [1]) def test_likelihoods_equal_priors(self): """likelihoods should equal Pr(D|H) if priors the same""" equal = [0.25, 0.25, 0.25,0.25] unequal = [0.5, 0.25, 0.125, 0.125] equal_answer = [1, 1, 1, 1] unequal_answer = [2, 1, 0.5, 0.5] for obs, exp in zip(likelihoods(equal, equal), equal_answer): self.assertFloatEqual(obs, exp) for obs, exp in zip(likelihoods(unequal, equal), unequal_answer): self.assertFloatEqual(obs, exp) def test_likelihoods_equal_evidence(self): """likelihoods should return vector of 1's if evidence equal for all""" equal = [0.25, 0.25, 0.25,0.25] unequal = [0.5, 0.25, 0.125, 0.125] equal_answer = [1, 1, 1, 1] unequal_answer = [2, 1, 0.5, 0.5] not_unity = [0.7, 0.7, 0.7, 0.7] for obs, exp in zip(likelihoods(equal, unequal), equal_answer): self.assertFloatEqual(obs, exp) #should be the same if evidences don't sum to 1 for obs, exp in zip(likelihoods(not_unity, unequal), equal_answer): self.assertFloatEqual(obs, exp) def test_likelihoods_unequal_evidence(self): """likelihoods should update based on weighted sum if evidence unequal""" not_unity = [1, 0.5, 0.25, 0.25] unequal = [0.5, 0.25, 0.125, 0.125] products = [1.4545455, 0.7272727, 0.3636364, 0.3636364] #if priors and evidence both unequal, likelihoods should change #(calculated using StarCalc) for obs, exp in zip(likelihoods(not_unity, unequal), products): self.assertFloatEqual(obs, exp) def test_posteriors_unequal_lists(self): """posteriors should raise ValueError if input lists unequal lengths""" self.assertRaises(ValueError, posteriors, [1, 2, 3], [1]) def test_posteriors_good_data(self): """posteriors should return products of paired list elements""" first = [0, 0.25, 0.5, 1, 0.25] second = [0.25, 0.5, 0, 0.1, 1] product = [0, 0.125, 0, 0.1, 0.25] for obs, exp in zip(posteriors(first, second), product): self.assertFloatEqual(obs, exp) class BayesUpdateTests(TestCase): """Tests implementation of Bayes calculations""" def setUp(self): first = [0.25, 0.25, 0.25] second = [0.1, 0.75, 0.3] third = [0.95, 1e-10, 0.2] fourth = [0.01, 0.9, 0.1] bad = [1, 2, 1, 1, 1] self.bad = [first, bad, second, third] self.test = [first, second, third, fourth] self.permuted = [fourth, first, third, second] self.deleted = [second, fourth, third] self.extra = [first, second, first, third, first, fourth, first] #BEWARE: low precision in second item, so need to adjust threshold #for assertFloatEqual accordingly (and use assertFloatEqualAbs). self.result = [0.136690646154, 0.000000009712, 0.863309344133] def test_bayes_updates_bad_data(self): """bayes_updates should raise ValueError on unequal-length lists""" self.assertRaises(ValueError, bayes_updates, self.bad) def test_bayes_updates_good_data(self): """bayes_updates should match hand calculations of probability updates""" #result for first -> fourth calculated by hand for obs, exp in zip(bayes_updates(self.test), self.result): self.assertFloatEqualAbs(obs, exp, 1e-11) def test_bayes_updates_permuted(self): """bayes_updates should not be affected by order of inputs""" for obs, exp in zip(bayes_updates(self.permuted), self.result): self.assertFloatEqualAbs(obs, exp, 1e-11) def test_bayes_update_nondiscriminating(self): """bayes_updates should be unaffected by extra nondiscriminating data""" #deletion of non-discriminating evidence should not affect result for obs, exp in zip(bayes_updates(self.deleted), self.result): self.assertFloatEqualAbs(obs, exp, 1e-11) #additional non-discriminating evidence should not affect result for obs, exp in zip(bayes_updates(self.extra), self.result): self.assertFloatEqualAbs(obs, exp, 1e-11) class StatTests(TestCase): """Tests that the t and z tests are implemented correctly""" def setUp(self): self.x = [ 7.33, 7.49, 7.27, 7.93, 7.56, 7.81, 7.46, 6.94, 7.49, 7.44, 7.95, 7.47, 7.04, 7.10, 7.64, ] self.y = [ 7.53, 7.70, 7.46, 8.21, 7.81, 8.01, 7.72, 7.13, 7.68, 7.66, 8.11, 7.66, 7.20, 7.25, 7.79, ] def test_t_paired_2tailed(self): """t_paired should match values from Sokal & Rohlf p 353""" x, y = self.x, self.y #check value of t and the probability for 2-tailed self.assertFloatEqual(t_paired(y, x)[0], 19.7203, 1e-4) self.assertFloatEqual(t_paired(y, x)[1], 1.301439e-11, 1e-4) def test_t_paired_no_variance(self): """t_paired should return None if lists are invariant""" x = [1, 1, 1] y = [0, 0, 0] self.assertEqual(t_paired(x,x), (None, None)) self.assertEqual(t_paired(x,y), (None, None)) def test_t_paired_1tailed(self): """t_paired should match pre-calculated 1-tailed values""" x, y = self.x, self.y #check probability for 1-tailed low and high self.assertFloatEqual( t_paired(y, x, "low")[1], 1-(1.301439e-11/2), 1e-4) self.assertFloatEqual( t_paired(x, y, "high")[1], 1-(1.301439e-11/2), 1e-4) self.assertFloatEqual( t_paired(y, x, "high")[1], 1.301439e-11/2, 1e-4) self.assertFloatEqual( t_paired(x, y, "low")[1], 1.301439e-11/2, 1e-4) def test_t_paired_specific_difference(self): """t_paired should allow a specific difference to be passed""" x, y = self.x, self.y #difference is 0.2, so test should be non-significant if 0.2 passed self.failIf(t_paired(y, x, exp_diff=0.2)[0] > 1e-10) #same, except that reversing list order reverses sign of difference self.failIf(t_paired(x, y, exp_diff=-0.2)[0] > 1e-10) #check that there's no significant difference from the true mean self.assertFloatEqual( t_paired(y, x,exp_diff=0.2)[1], 1, 1e-4) def test_t_paired_bad_data(self): """t_paired should raise ValueError on lists of different lengths""" self.assertRaises(ValueError, t_paired, self.y, [1, 2, 3]) def test_t_two_sample(self): """t_two_sample should match example on p.225 of Sokal and Rohlf""" I = array([7.2, 7.1, 9.1, 7.2, 7.3, 7.2, 7.5]) II = array([8.8, 7.5, 7.7, 7.6, 7.4, 6.7, 7.2]) self.assertFloatEqual(t_two_sample(I, II), (-0.1184, 0.45385 * 2), 0.001) def test_t_two_sample_no_variance(self): """t_two_sample should return None if lists are invariant""" x = array([1, 1, 1]) y = array([0, 0, 0]) self.assertEqual(t_two_sample(x,x), (None, None)) self.assertEqual(t_two_sample(x,y), (None, None)) def test_t_one_sample(self): """t_one_sample results should match those from R""" x = array(range(-5,5)) y = array(range(-1,10)) self.assertFloatEqualAbs(t_one_sample(x), (-0.5222, 0.6141), 1e-4) self.assertFloatEqualAbs(t_one_sample(y), (4, 0.002518), 1e-4) #do some one-tailed tests as well self.assertFloatEqualAbs(t_one_sample(y, tails='low'),(4, 0.9987),1e-4) self.assertFloatEqualAbs(t_one_sample(y,tails='high'),(4,0.001259),1e-4) def test_t_two_sample_switch(self): """t_two_sample should call t_one_observation if 1 item in sample.""" sample = array([4.02, 3.88, 3.34, 3.87, 3.18]) x = array([3.02]) self.assertFloatEqual(t_two_sample(x,sample),(-1.5637254,0.1929248)) self.assertFloatEqual(t_two_sample(sample, x),(-1.5637254,0.1929248)) #can't do the test if both samples have single item self.assertEqual(t_two_sample(x,x), (None, None)) def test_t_one_observation(self): """t_one_observation should match p. 228 of Sokal and Rohlf""" sample = array([4.02, 3.88, 3.34, 3.87, 3.18]) x = 3.02 #note that this differs after the 3rd decimal place from what's in the #book, because Sokal and Rohlf round their intermediate steps... self.assertFloatEqual(t_one_observation(x,sample),\ (-1.5637254,0.1929248)) def test_reverse_tails(self): """reverse_tails should return 'high' if tails was 'low' or vice versa""" self.assertEqual(reverse_tails('high'), 'low') self.assertEqual(reverse_tails('low'), 'high') self.assertEqual(reverse_tails(None), None) self.assertEqual(reverse_tails(3), 3) def test_tail(self): """tail should return prob/2 if test is true, or 1-(prob/2) if false""" self.assertFloatEqual(tail(0.25, True), 0.125) self.assertFloatEqual(tail(0.25, False), 0.875) self.assertFloatEqual(tail(1, True), 0.5) self.assertFloatEqual(tail(1, False), 0.5) self.assertFloatEqual(tail(0, True), 0) self.assertFloatEqual(tail(0, False), 1) def test_z_test(self): """z_test should give correct values""" sample = array([1,2,3,4,5]) self.assertFloatEqual(z_test(sample, 3, 1), (0,1)) self.assertFloatEqual(z_test(sample, 3, 2, 'high'), (0,0.5)) self.assertFloatEqual(z_test(sample, 3, 2, 'low'), (0,0.5)) #check that population mean and variance, and tails, can be set OK. self.assertFloatEqual(z_test(sample, 0, 1), (6.7082039324993694, \ 1.9703444711798951e-11)) self.assertFloatEqual(z_test(sample, 1, 10), (0.44721359549995793, \ 0.65472084601857694)) self.assertFloatEqual(z_test(sample, 1, 10, 'high'), \ (0.44721359549995793, 0.65472084601857694/2)) self.assertFloatEqual(z_test(sample, 1, 10, 'low'), \ (0.44721359549995793, 1-(0.65472084601857694/2))) class CorrelationTests(TestCase): """Tests of correlation coefficients.""" def test_correlation(self): """Correlations and significance should match R's cor.test()""" x = [1,2,3,5] y = [0,0,0,0] z = [1,1,1,1] a = [2,4,6,8] b = [1.5, 1.4, 1.2, 1.1] c = [15, 10, 5, 20] bad = [1,2,3] #originally gave r = 1.0000000002 self.assertFloatEqual(correlation(x,x), (1, 0)) self.assertFloatEqual(correlation(x,y), (0,1)) self.assertFloatEqual(correlation(y,z), (0,1)) self.assertFloatEqualAbs(correlation(x,a), (0.9827076, 0.01729), 1e-5) self.assertFloatEqualAbs(correlation(x,b), (-0.9621405, 0.03786), 1e-5) self.assertFloatEqualAbs(correlation(x,c), (0.3779645, 0.622), 1e-3) self.assertEqual(correlation(bad,bad), (1, 0)) def test_correlation_matrix(self): """Correlations in matrix should match values from R""" a = [2,4,6,8] b = [1.5, 1.4, 1.2, 1.1] c = [15, 10, 5, 20] m = correlation_matrix([a,b,c]) self.assertFloatEqual(m[0,0], [1.0]) self.assertFloatEqual([m[1,0], m[1,1]], [correlation(b,a)[0], 1.0]) self.assertFloatEqual(m[2], [correlation(c,a)[0], correlation(c,b)[0], \ 1.0]) class Ftest(TestCase): """Tests for the F test""" def test_f_value(self): """f_value: should calculate the correct F value if possible""" a = array([1,3,5,7,9,8,6,4,2]) b = array([5,4,6,3,7,6,4,5]) self.assertEqual(f_value(a,b), (8,7,4.375)) self.assertFloatEqual(f_value(b,a), (7,8,0.2285714)) too_short = array([4]) self.assertRaises(ValueError, f_value, too_short, b) def test_f_two_sample(self): """f_two_sample should match values from R""" #The expected values in this test are obtained through R. #In R the F test is var.test(x,y) different alternative hypotheses #can be specified (two sided, less, or greater). #The vectors are random samples from a particular normal distribution #(mean and sd specified). #a: 50 elem, mean=0 sd=1 a = [-0.70701689, -1.24788845, -1.65516470, 0.10443876, -0.48526915, -0.71820656, -1.02603596, 0.03975982, -2.23404324, -0.21509363, 0.08438468, -0.01970062, -0.67907971, -0.89853667, 1.11137131, 0.05960496, -1.51172084, -0.79733957, -1.60040659, 0.80530639, -0.81715836, -0.69233474, 0.95750665, 0.99576429, -1.61340216, -0.43572590, -1.50862327, 0.92847551, -0.68382338, -1.12523522, -0.09147488, 0.66756023, -0.87277588, -1.36539039, -0.11748707, -1.63632578, -0.31343078, -0.28176086, 0.33854483, -0.51785630, 2.25360559, -0.80761191, 1.18983499, 0.57080342, -1.44601700, -0.53906955, -0.01975266, -1.37147915, -0.31537616, 0.26877544] #b: 50 elem, mean=0, sd=1.2 b=[0.081418743, 0.276571612, -1.864316504, 0.675213612, -0.769202643, 0.140372825, -1.426250184, 0.058617884, -0.819287409, -0.007701916, -0.782722020, -0.285891593, 0.661980419, 0.383225191, 0.622444946, -0.192446150, 0.297150571, 0.408896059, -0.167359383, -0.552381362, 0.982168338, 1.439730446, 1.967616101, -0.579607307, 1.095590943, 0.240591302, -1.566937143, -0.199091349, -1.232983905, 0.362378169, 1.166061081, -0.604676222, -0.536560206, -0.303117595, 1.519222792, -0.319146503, 2.206220810, -0.566351124, -0.720397392, -0.452001377, 0.250890097, 0.320685395, -1.014632725, -3.010346273, -1.703955054, 0.592587381, -1.237451255, 0.172243366, -0.452641122, -0.982148581] #c: 60 elem, mean=5, sd=1 c=[4.654329, 5.242129, 6.272640, 5.781779, 4.391241, 3.800752, 4.559463, 4.318922, 3.243020, 5.121280, 4.126385, 5.541131, 4.777480, 5.646913, 6.972584, 3.817172, 6.128700, 4.731467, 6.762068, 5.082983, 5.298511, 5.491125, 4.532369, 4.265552, 5.697317, 5.509730, 2.935704, 4.507456, 3.786794, 5.548383, 3.674487, 5.536556, 5.297847, 2.439642, 4.759836, 5.114649, 5.986774, 4.517485, 4.579208, 4.579374, 2.502890, 5.190955, 5.983194, 6.766645, 4.905079, 4.214273, 3.950364, 6.262393, 8.122084, 6.330007, 4.767943, 5.194029, 3.503136, 6.039079, 4.485647, 6.116235, 6.302268, 3.596693, 5.743316, 6.860152] #d: 30 elem, mean=0, sd =0.05 d=[ 0.104517366, 0.023039678, 0.005579091, 0.052928250, 0.020724823, -0.060823243, -0.019000890, -0.064133996, -0.016321594, -0.008898334, -0.027626992, -0.051946186, 0.085269587, -0.031190678, 0.065172938, -0.054628573, 0.019257306, -0.032427056, -0.058767356, 0.030927400, 0.052247357, -0.042954937, 0.031842104, 0.094130522, -0.024828465, 0.011320453, -0.016195062, 0.015631245, -0.050335598, -0.031658335] a,b,c,d = map(array,[a,b,c,d]) self.assertEqual(map(len,[a,b,c,d]), [50, 50, 60, 30]) #allowed error. This big, because results from R #are rounded at 4 decimals error = 1e-4 self.assertFloatEqual(f_two_sample(a,a), (49, 49, 1, 1), eps=error) self.assertFloatEqual(f_two_sample(a,b), (49, 49, 0.8575, 0.5925), eps=error) self.assertFloatEqual(f_two_sample(b,a), (49, 49, 1.1662, 0.5925), eps=error) self.assertFloatEqual(f_two_sample(a,b, tails='low'), (49, 49, 0.8575, 0.2963), eps=error) self.assertFloatEqual(f_two_sample(a,b, tails='high'), (49, 49, 0.8575, 0.7037), eps=error) self.assertFloatEqual(f_two_sample(a,c), (49, 59, 0.6587, 0.1345), eps=error) #p value very small, so first check df's and F value self.assertFloatEqualAbs(f_two_sample(d,a, tails='low')[0:3], (29, 49, 0.0028), eps=error) assert f_two_sample(d,a, tails='low')[3] < 2.2e-16 #p value def test_MonteCarloP(self): """MonteCarloP calcs a p-value from a val and list of random vals""" val = 3.0 random_vals = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0] #test for "high" tail (larger values than expected by chance) p_val = MonteCarloP(val, random_vals, 'high') self.assertEqual(p_val, 0.7) #test for "low" tail (smaller values than expected by chance) p_val = MonteCarloP(val, random_vals, 'low') self.assertEqual(p_val, 0.4) def test_permute_2d(self): """permute_2d permutes rows and cols of a matrix.""" a = reshape(arange(9), (3,3)) self.assertEqual(permute_2d(a, [0,1,2]), a) self.assertEqual(permute_2d(a, [2,1,0]), \ array([[8,7,6],[5,4,3],[2,1,0]])) self.assertEqual(permute_2d(a, [1,2,0]), \ array([[4,5,3],[7,8,6],[1,2,0]])) def test_mantel(self): """mantel should be significant for same matrix, not for random""" a = reshape(arange(25), (5,5)) b = a.copy() b[-1,-1] = 26 #slight change m = mantel(a, b, 1000) #closely related -- should be significant assert m < 0.05 c = reshape(ones(25), (5,5)) c[-1,-1] = 3 #not related -- should not be significant m = mantel(a,c,1000) assert m > 0.3, m class KendallTests(TestCase): """check accuracy of Kendall tests against values from R""" def do_test(self, x, y, alt_expecteds): """conducts the tests for each alternate hypothesis against expecteds""" for alt, exp_p, exp_tau in alt_expecteds: tau, p_val = kendall_correlation(x, y, alt=alt, warn=False) self.assertFloatEqual(tau, exp_tau, eps=1e-3) self.assertFloatEqual(p_val, exp_p, eps=1e-3) def test_exact_calcs(self): """calculations of exact probabilities should match R""" x = (44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1) y = ( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8) expecteds = [["gt", 0.05972, 0.4444444], ["lt", 0.9624, 0.4444444], ["ts", 0.1194, 0.4444444]] self.do_test(x,y,expecteds) def test_with_ties(self): """tied values calculated from normal approx""" # R example with ties in x x = (44.4, 45.9, 41.9, 53.3, 44.4, 44.1, 50.7, 45.2, 60.1) y = ( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8) expecteds = [#["gt", 0.05793, 0.4225771], ["lt", 0.942, 0.4225771], ["ts", 0.1159, 0.4225771]] self.do_test(x,y,expecteds) # R example with ties in y x = (44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1) y = ( 2.6, 3.1, 2.5, 5.0, 3.1, 4.0, 5.2, 2.8, 3.8) expecteds = [["gt", 0.03737, 0.4789207], ["lt", 0.9626, 0.4789207], ["ts", 0.07474, 0.4789207]] self.do_test(x,y,expecteds) # R example with ties in x and y x = (44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 44.4, 60.1) y = ( 2.6, 3.6, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8) expecteds=[["gt", 0.02891, 0.5142857], ["lt", 0.971, 0.5142857], ["ts", 0.05782, 0.5142857]] self.do_test(x,y,expecteds) def test_bigger_vectors(self): """docstring for test_bigger_vectors""" # q < expansion x= (0.118583104633, 0.227860069338, 0.143856130991, 0.935362617582, 0.0471303856799, 0.659819202174, 0.739247965907, 0.268929000278, 0.848250568194, 0.307764819102, 0.733949480141, 0.271662210481, 0.155903098872) y= (0.749762144455, 0.407571703468, 0.934176427266, 0.188638794706, 0.184844781493, 0.391485553856, 0.735504815302, 0.363655952442, 0.18489971978, 0.851075466765, 0.139932273818, 0.333675110224, 0.570250937033) expecteds = [["gt", 0.9183, -0.2820513], ["lt", 0.1022, -0.2820513], ["ts", 0.2044, -0.2820513]] self.do_test(x,y,expecteds) # q > expansion x= (0.2602556958, 0.441506392849, 0.930624643531, 0.728461775775, 0.234341774892, 0.725677256368, 0.354788882728, 0.475882541956, 0.347533553428, 0.608578046857, 0.144697962102, 0.784502692164, 0.872607603407) y= (0.753056395718, 0.454332072011, 0.791882395707, 0.622853579015, 0.127030232518, 0.232086215578, 0.586604349918, 0.0139051260749, 0.579079370051, 0.0550643809812, 0.94798878249, 0.318410679439, 0.86725134615) expecteds = [["gt", 0.4762, 0.02564103], ["lt", 0.5711, 0.02564103], ["ts", 0.9524, 0.02564103]] self.do_test(x,y,expecteds) class TestDistMatrixPermutationTest(TestCase): """Tests of distance_matrix_permutation_test""" def setUp(self): """sets up variables for testing""" self.matrix = array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]]) self.cells = [(0,1), (1,3)] self.cells2 = [(0,2), (2,3)] def test_get_ltm_cells(self): "get_ltm_cells converts indices to be below the diagonal" cells = [(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)] result = get_ltm_cells(cells) self.assertEqual(result, [(2, 0), (1, 0), (2, 1)]) cells = [(0,1),(0,2)] result = get_ltm_cells(cells) self.assertEqual(result, [(2, 0), (1, 0)]) def test_get_values_from_matrix(self): """get_values_from_matrix returns the special and other values from matrix""" matrix = self.matrix cells = [(1,0),(0,1),(2,0),(2,1)] #test that works for a symmetric matrix cells_sym = get_ltm_cells(cells) special_vals, other_vals = get_values_from_matrix(matrix, cells_sym,\ cells2=None, is_symmetric=True) special_vals.sort() other_vals.sort() self.assertEqual(special_vals, [5,9,10]) self.assertEqual(other_vals, [13,14,15]) #test that work for a non symmetric matrix special_vals, other_vals = get_values_from_matrix(matrix, cells,\ cells2=None, is_symmetric=False) special_vals.sort() other_vals.sort() self.assertEqual(special_vals, [2,5,9,10]) self.assertEqual(other_vals, [1,3,4,6,7,8,11,12,13,14,15,16]) #test that works on a symmetric matrix when cells2 is defined cells2 = [(3,0),(3,2),(0,3)] cells2_sym = get_ltm_cells(cells2) special_vals, other_vals = get_values_from_matrix(matrix, cells_sym,\ cells2=cells2_sym, is_symmetric=True) special_vals.sort() other_vals.sort() self.assertEqual(special_vals, [5,9,10]) self.assertEqual(other_vals, [13,15]) #test that works when cells2 is defined and not symmetric special_vals, other_vals = get_values_from_matrix(matrix, cells, cells2=cells2,\ is_symmetric=False) special_vals.sort() other_vals.sort() self.assertEqual(special_vals, [2,5,9,10]) self.assertEqual(other_vals, [4,13,15]) def test_distance_matrix_permutation_test_non_symmetric(self): """ evaluate empirical p-values for a non symmetric matrix To test the empirical p-values, we look at a simple 3x3 matrix b/c it is easy to see what t score every permutation will generate -- there's only 6 permutations. Running dist_matrix_test with n=1000, we expect that each permutation will show up 160 times, so we know how many times to expect to see more extreme t scores. We therefore know what the empirical p-values will be. (n=1000 was chosen empirically -- smaller values seem to lead to much more frequent random failures.) """ def make_result_list(*args, **kwargs): return [distance_matrix_permutation_test(*args,**kwargs)[2] \ for i in range(10)] m = arange(9).reshape((3,3)) n = 100 # looks at each possible permutation n times -- # compare first row to rest r = make_result_list(m, [(0,0),(0,1),(0,2)],n=n,is_symmetric=False) self.assertSimilarMeans(r, 0./6.) r = make_result_list(m, [(0,0),(0,1),(0,2)],n=n,is_symmetric=False,\ tails='high') self.assertSimilarMeans(r, 4./6.) r = make_result_list(m, [(0,0),(0,1),(0,2)],n=n,is_symmetric=False,\ tails='low') self.assertSimilarMeans(r, 0./6.) # looks at each possible permutation n times -- # compare last row to rest r = make_result_list(m, [(2,0),(2,1),(2,2)],n=n,is_symmetric=False) self.assertSimilarMeans(r, 0./6.) r = make_result_list(m, [(2,0),(2,1),(2,2)],n=n,is_symmetric=False,\ tails='high') self.assertSimilarMeans(r, 0./6.) r = make_result_list(m, [(2,0),(2,1),(2,2)],n=n,is_symmetric=False,\ tails='low') self.assertSimilarMeans(r, 4./6.) def test_distance_matrix_permutation_test_symmetric(self): """ evaluate empirical p-values for symmetric matrix See test_distance_matrix_permutation_test_non_symmetric doc string for a description of how this test works. """ def make_result_list(*args, **kwargs): return [distance_matrix_permutation_test(*args)[2] for i in range(10)] m = array([[0,1,3],[1,2,4],[3,4,5]]) # looks at each possible permutation n times -- # compare first row to rest n = 100 # looks at each possible permutation n times -- # compare first row to rest r = make_result_list(m, [(0,0),(0,1),(0,2)],n=n) self.assertSimilarMeans(r, 0./6.) r = make_result_list(m, [(0,0),(0,1),(0,2)],n=n,tails='high') self.assertSimilarMeans(r, 0.77281447417149496,0) r = make_result_list(m, [(0,0),(0,1),(0,2)],n=n,tails='low') self.assertSimilarMeans(r, 4./6.) ## The following lines are not part of the test code, but are useful in ## figuring out what t-scores all of the permutations will yield. #permutes = [[0, 1, 2], [0, 2, 1], [1, 0, 2],\ # [1, 2, 0], [2, 0, 1], [2, 1, 0]] #results = [] #for p in permutes: # p_m = permute_2d(m,p) # results.append(t_two_sample(\ # [p_m[0,1],p_m[0,2]],[p_m[2,1]],tails='high')) #print results def test_distance_matrix_permutation_test_alt_stat(self): def fake_stat_test(a,b,tails=None): return 42.,42. m = array([[0,1,3],[1,2,4],[3,4,5]]) self.assertEqual(distance_matrix_permutation_test(m,\ [(0,0),(0,1),(0,2)],n=5,f=fake_stat_test),(42.,42.,0.)) def test_distance_matrix_permutation_test_return_scores(self): """ return_scores=True functions as expected """ # use alt statistical test to make results simple def fake_stat_test(a,b,tails=None): return 42.,42. m = array([[0,1,3],[1,2,4],[3,4,5]]) self.assertEqual(distance_matrix_permutation_test(\ m,[(0,0),(0,1),(0,2)],\ n=5,f=fake_stat_test,return_scores=True),(42.,42.,0.,[42.]*5)) def test_ANOVA_one_way(self): """ANOVA one way returns same values as ANOVA on a stats package """ g1 = Numbers([10.0, 11.0, 10.0, 5.0, 6.0]) g2 = Numbers([1.0, 2.0, 3.0, 4.0, 1.0, 2.0]) g3 = Numbers([6.0, 7.0, 5.0, 6.0, 7.0]) i = [g1, g2, g3] dfn, dfd, F, between_MS, within_MS, group_means, prob = ANOVA_one_way(i) self.assertEqual(dfn, 2) self.assertEqual(dfd, 13) self.assertFloatEqual(F, 18.565450643776831) self.assertFloatEqual(between_MS, 55.458333333333343) self.assertFloatEqual(within_MS, 2.9871794871794868) self.assertFloatEqual(group_means, [8.4000000000000004, 2.1666666666666665, 6.2000000000000002]) self.assertFloatEqual(prob, 0.00015486238993089464) #execute tests if called from command line if __name__ == '__main__': main()