mvpa.tests.test

1 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 2 # vi: set ft=python sts=4 ts=4 sw=4 et: 3 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 4 # 5 # See COPYING file distributed along with the PyMVPA package for the 6 # copyright and license terms. 7 # 8 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 9 """Unit tests for PyMVPA stats helpers""" 10 11 from mvpa.base import externals 12 from mvpa.clfs.stats import MCNullDist, FixedNullDist, NullDist 13 from mvpa.datasets import Dataset 14 from mvpa.measures.glm import GLM 15 from mvpa.measures.anova import OneWayAnova, CompoundOneWayAnova 16 from mvpa.misc.fx import doubleGammaHRF, singleGammaHRF 17 from tests_warehouse import * 18 from mvpa import cfg 19 20 # Prepare few distributions to test 21 #kwargs = {'permutations':10, 'tail':'any'} 22 nulldist_sweep = [ MCNullDist(permutations=10, tail='any'), 23 MCNullDist(permutations=10, tail='right')] 24 if externals.exists('scipy'): 25 from mvpa.support.stats import scipy 26 nulldist_sweep += [ MCNullDist(scipy.stats.norm, permutations=10, 27 tail='any'), 28 MCNullDist(scipy.stats.norm, permutations=10, 29 tail='right'), 30 MCNullDist(scipy.stats.expon, permutations=10, 31 tail='right'), 32 FixedNullDist(scipy.stats.norm(0, 0.1), tail='any'), 33 FixedNullDist(scipy.stats.norm(0, 0.1), tail='right'), 34 scipy.stats.norm(0, 0.1) 35 ]

36 37 -class StatsTests(unittest.TestCase):

38 """Unittests for various statistics""" 39

40 - def testChiSquare(self):

41 """Test chi-square distribution""" 42 if not externals.exists('scipy'): 43 return 44 45 from mvpa.misc.stats import chisquare 46 47 # test equal distribution 48 tbl = N.array([[5, 5], [5, 5]]) 49 chi, p = chisquare(tbl) 50 self.failUnless( chi == 0.0 ) 51 self.failUnless( p == 1.0 ) 52 53 # test non-equal distribution 54 tbl = N.array([[4, 0], [0, 4]]) 55 chi, p = chisquare(tbl) 56 self.failUnless(chi == 8.0) 57 self.failUnless(p < 0.05)

58 59 60 @sweepargs(null=nulldist_sweep[1:])

61 - def testNullDistProb(self, null):

62 """Testing null dist probability""" 63 if not isinstance(null, NullDist): 64 return 65 ds = datasets['uni2small'] 66 67 null.fit(OneWayAnova(), ds) 68 69 # check reasonable output. 70 # p-values for non-bogus features should significantly different, 71 # while bogus (0) not 72 prob = null.p([3, 0, 0, 0, 0, N.nan]) 73 self.failUnless(N.abs(prob[0]) < 0.01) 74 if cfg.getboolean('tests', 'labile', default='yes'): 75 self.failUnless((N.abs(prob[1:]) > 0.05).all(), 76 msg="Bogus features should have insignificant p." 77 " Got %s" % (N.abs(prob[1:]),)) 78 79 # has to have matching shape 80 if not isinstance(null, FixedNullDist): 81 # Fixed dist is univariate ATM so it doesn't care 82 # about dimensionality and gives 1 output value 83 self.failUnlessRaises(ValueError, null.p, [5, 3, 4])

84 85

86 - def testNullDistProbAny(self):

87 """Test 'any' tail statistics estimation""" 88 if not externals.exists('scipy'): 89 return 90 91 # test 'any' mode 92 from mvpa.measures.corrcoef import CorrCoef 93 ds = datasets['uni2small'] 94 95 null = MCNullDist(permutations=10, tail='any') 96 null.fit(CorrCoef(), ds) 97 98 # 100 and -100 should both have zero probability on their respective 99 # tails 100 self.failUnless(null.p([-100, 0, 0, 0, 0, 0])[0] == 0) 101 self.failUnless(null.p([100, 0, 0, 0, 0, 0])[0] == 0) 102 103 # same test with just scalar measure/feature 104 null.fit(CorrCoef(), ds.selectFeatures([0])) 105 self.failUnless(null.p(-100) == 0) 106 self.failUnless(null.p(100) == 0)

107 108 109 @sweepargs(nd=nulldist_sweep)

110 - def testDatasetMeasureProb(self, nd):

111 """Test estimation of measures statistics""" 112 if not externals.exists('scipy'): 113 # due to null_t requirement 114 return 115 116 ds = datasets['uni2medium'] 117 118 m = OneWayAnova(null_dist=nd, enable_states=['null_t']) 119 score = m(ds) 120 121 score_nonbogus = N.mean(score[ds.nonbogus_features]) 122 score_bogus = N.mean(score[ds.bogus_features]) 123 # plausability check 124 self.failUnless(score_bogus < score_nonbogus) 125 126 null_prob_nonbogus = m.null_prob[ds.nonbogus_features] 127 null_prob_bogus = m.null_prob[ds.bogus_features] 128 129 self.failUnless((null_prob_nonbogus < 0.05).all(), 130 msg="Nonbogus features should have a very unlikely value. Got %s" 131 % null_prob_nonbogus) 132 133 # the others should be a lot larger 134 self.failUnless(N.mean(N.abs(null_prob_bogus)) > 135 N.mean(N.abs(null_prob_nonbogus))) 136 137 if isinstance(nd, MCNullDist): 138 # MCs are not stable with just 10 samples... so lets skip them 139 return 140 141 if cfg.getboolean('tests', 'labile', default='yes'): 142 # Failed on c94ec26eb593687f25d8c27e5cfdc5917e352a69 143 # with MVPA_SEED=833393575 144 self.failUnless((N.abs(m.null_t[ds.nonbogus_features]) >= 5).all(), 145 msg="Nonbogus features should have high t-score. Got %s" 146 % (m.null_t[ds.nonbogus_features])) 147 148 self.failUnless((N.abs(m.null_t[ds.bogus_features]) < 4).all(), 149 msg="Bogus features should have low t-score. Got (t,p,sens):%s" 150 % (zip(m.null_t, m.null_prob, score)))

151 152

153 - def testNegativeT(self):

154 """Basic testing of the sign in p and t scores 155 """ 156 from mvpa.measures.base import FeaturewiseDatasetMeasure 157 158 class BogusMeasure(FeaturewiseDatasetMeasure): 159 """Just put high positive into first 2 features, and high 160 negative into 2nd two 161 """ 162 def _call(self, dataset): 163 """just a little helper... pylint shut up!""" 164 res = N.random.normal(size=(dataset.nfeatures,)) 165 res[0] = res[1] = 100 166 res[2] = res[3] = -100 167 return res

168 if not externals.exists('scipy'): 169 return 170 nd = FixedNullDist(scipy.stats.norm(0, 0.1), tail='any') 171 m = BogusMeasure(null_dist=nd, enable_states=['null_t']) 172 ds = datasets['uni2small'] 173 score = m(ds) 174 t, p = m.null_t, m.null_prob 175 self.failUnless((p>=0).all()) 176 self.failUnless((t[:2] > 0).all()) 177 self.failUnless((t[2:4] < 0).all()) 178 179

180 - def testMatchDistribution(self):

181 """Some really basic testing for matchDistribution 182 """ 183 if not externals.exists('scipy'): 184 return 185 from mvpa.clfs.stats import matchDistribution, rv_semifrozen 186 187 data = datasets['uni2small'].samples[:, 1] 188 189 # Lets test ad-hoc rv_semifrozen 190 floc = rv_semifrozen(scipy.stats.norm, loc=0).fit(data) 191 self.failUnless(floc[0] == 0) 192 193 fscale = rv_semifrozen(scipy.stats.norm, scale=1.0).fit(data) 194 self.failUnless(fscale[1] == 1) 195 196 flocscale = rv_semifrozen(scipy.stats.norm, loc=0, scale=1.0).fit(data) 197 self.failUnless(flocscale[1] == 1 and flocscale[0] == 0) 198 199 full = scipy.stats.norm.fit(data) 200 for res in [floc, fscale, flocscale, full]: 201 self.failUnless(len(res) == 2) 202 203 data_mean = N.mean(data) 204 for loc in [None, data_mean]: 205 for test in ['p-roc', 'kstest']: 206 # some really basic testing 207 matched = matchDistribution( 208 data=data, 209 distributions = ['scipy', 210 ('norm', 211 {'name': 'norm_fixed', 212 'loc': 0.2, 213 'scale': 0.3})], 214 test=test, loc=loc, p=0.05) 215 # at least norm should be in there 216 names = [m[2] for m in matched] 217 if test == 'p-roc': 218 if cfg.getboolean('tests', 'labile', default='yes'): 219 # we can guarantee that only for norm_fixed 220 self.failUnless('norm' in names) 221 self.failUnless('norm_fixed' in names) 222 inorm = names.index('norm_fixed') 223 # and it should be at least in the first 224 # 30 best matching ;-) 225 self.failUnless(inorm <= 30)

226 227

228 - def testRDistStability(self):

229 """Test either rdist distribution performs nicely 230 """ 231 if not externals.exists('scipy'): 232 return 233 234 try: 235 # actually I haven't managed to cause this error 236 scipy.stats.rdist(1.32, 0, 1).pdf(-1.0+N.finfo(float).eps) 237 except Exception, e: 238 self.fail('Failed to compute rdist.pdf due to numeric' 239 ' loss of precision. Exception was %s' % e) 240 241 try: 242 # this one should fail with recent scipy with error 243 # ZeroDivisionError: 0.0 cannot be raised to a negative power 244 245 # XXX: There is 1 more bug in etch's scipy.stats or numpy 246 # (vectorize), so I have to put 2 elements in the 247 # queried x's, otherwise it 248 # would puke. But for now that fix is not here 249 # 250 # value = scipy.stats.rdist(1.32, 0, 1).cdf( 251 # [-1.0+N.finfo(float).eps, 0]) 252 # 253 # to cause it now just run this unittest only with 254 # nosetests -s test_stats:StatsTests.testRDistStability 255 256 # NB: very cool way to store the trace of the execution 257 #import pydb 258 #pydb.debugger(dbg_cmds=['bt', 'l', 's']*300 + ['c']) 259 scipy.stats.rdist(1.32, 0, 1).cdf(-1.0+N.finfo(float).eps) 260 except IndexError, e: 261 self.fail('Failed due to bug which leads to InvalidIndex if only' 262 ' scalar is provided to cdf') 263 except Exception, e: 264 self.fail('Failed to compute rdist.cdf due to numeric' 265 ' loss of precision. Exception was %s' % e) 266 267 v = scipy.stats.rdist(10000, 0, 1).cdf([-0.1]) 268 self.failUnless(v>=0, v<=1)

269 270

271 - def testAnova(self):

272 """Do some extended testing of OneWayAnova 273 274 in particular -- compound estimation 275 """ 276 277 m = OneWayAnova() # default must be not compound ? 278 mc = CompoundOneWayAnova(combiner=None) 279 ds = datasets['uni2medium'] 280 281 # For 2 labels it must be identical for both and equal to 282 # simple OneWayAnova 283 a, ac = m(ds), mc(ds) 284 285 self.failUnless(a.shape == (ds.nfeatures,)) 286 self.failUnless(ac.shape == (ds.nfeatures, len(ds.uniquelabels))) 287 288 self.failUnless((ac[:, 0] == ac[:, 1]).all()) 289 self.failUnless((a == ac[:, 1]).all()) 290 291 ds = datasets['uni4large'] 292 ac = mc(ds) 293 if cfg.getboolean('tests', 'labile', default='yes'): 294 # All non-bogus features must be high for a corresponding feature 295 self.failUnless((ac[(N.array(ds.nonbogus_features), 296 N.arange(4))] >= 1).all())

297 298

299 - def testGLM(self):

300 """Test GLM 301 """ 302 # play fmri 303 # full-blown HRF with initial dip and undershoot ;-) 304 hrf_x = N.linspace(0, 25, 250) 305 hrf = doubleGammaHRF(hrf_x) - singleGammaHRF(hrf_x, 0.8, 1, 0.05) 306 307 # come up with an experimental design 308 samples = 1800 309 fast_er_onsets = N.array([10, 200, 250, 500, 600, 900, 920, 1400]) 310 fast_er = N.zeros(samples) 311 fast_er[fast_er_onsets] = 1 312 313 # high resolution model of the convolved regressor 314 model_hr = N.convolve(fast_er, hrf)[:samples] 315 316 if not externals.exists('scipy'): 317 return 318 319 from scipy import signal 320 # downsample the regressor to fMRI resolution 321 tr = 2.0 322 model_lr = signal.resample(model_hr, 323 int(samples / tr / 10), 324 window='ham') 325 326 # generate artifical fMRI data: two voxels one is noise, one has 327 # something 328 baseline = 800.0 329 wsignal = baseline + 2 * model_lr + \ 330 N.random.randn(int(samples / tr / 10)) * 0.2 331 nsignal = baseline + N.random.randn(int(samples / tr / 10)) * 0.5 332 333 # build design matrix: bold-regressor and constant 334 X = N.array([model_lr, N.repeat(1, len(model_lr))]).T 335 336 # two 'voxel' dataset 337 data = Dataset(samples=N.array((wsignal, nsignal, nsignal)).T, labels=1) 338 339 # check GLM betas 340 glm = GLM(X, combiner=None) 341 betas = glm(data) 342 343 # betas for each feature and each regressor 344 self.failUnless(betas.shape == (data.nfeatures, X.shape[1])) 345 346 self.failUnless(N.absolute(betas[:, 1] - baseline < 10).all(), 347 msg="baseline betas should be huge and around 800") 348 349 self.failUnless(betas[0][0] > betas[1, 0], 350 msg="feature (with signal) beta should be larger than for noise") 351 352 if cfg.getboolean('tests', 'labile', default='yes'): 353 self.failUnless(N.absolute(betas[1, 0]) < 0.5) 354 self.failUnless(N.absolute(betas[0, 0]) > 1.0) 355 356 357 # check GLM zscores 358 glm = GLM(X, voi='zstat', combiner=None) 359 zstats = glm(data) 360 361 self.failUnless(zstats.shape == betas.shape) 362 363 self.failUnless((zstats[:, 1] > 1000).all(), 364 msg='constant zstats should be huge') 365 366 if cfg.getboolean('tests', 'labile', default='yes'): 367 self.failUnless(N.absolute(betas[0, 0]) > betas[1][0], 368 msg='with signal should have higher zstats')

369

370 371 -def suite():

372 """Create the suite""" 373 return unittest.makeSuite(StatsTests)

374 375 376 if __name__ == '__main__': 377 import runner 378

Source Code for Module mvpa.tests.test_stats