mvpa.clfs.stats

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 """Estimator for classifier error distributions.""" 10 11 __docformat__ = 'restructuredtext' 12 13 import numpy as N 14 15 from mvpa.base import externals, warning 16 from mvpa.misc.state import ClassWithCollections, StateVariable 17 18 if __debug__: 19 from mvpa.base import debug

20 21 22 -class Nonparametric(object):

23 """Non-parametric 1d distribution -- derives cdf based on stored values. 24 25 Introduced to complement parametric distributions present in scipy.stats. 26 """ 27

28 - def __init__(self, dist_samples, correction='clip'):

29 """ 30 :Parameters: 31 dist_samples : ndarray 32 Samples to be used to assess the distribution. 33 correction : {'clip'} or None, optional 34 Determines the behavior when .cdf is queried. If None -- no 35 correction is made. If 'clip' -- values are clipped to lie 36 in the range [1/(N+2), (N+1)/(N+2)] (simply because 37 non-parametric assessment lacks the power to resolve with 38 higher precision in the tails, so 'imagery' samples are 39 placed in each of the two tails). 40 """ 41 self._dist_samples = N.ravel(dist_samples) 42 self._correction = correction

43

44 - def __repr__(self):

45 return '%s(%r%s)' % ( 46 self.__class__.__name__, 47 self._dist_samples, 48 ('', ', correction=%r' % self._correction) 49 [int(self._correction != 'clip')])

50 51 @staticmethod

52 - def fit(dist_samples):

53 return [dist_samples]

54 55

56 - def cdf(self, x):

57 """Returns the cdf value at `x`. 58 """ 59 dist_samples = self._dist_samples 60 res = N.vectorize(lambda v:(dist_samples <= v).mean())(x) 61 if self._correction == 'clip': 62 nsamples = len(dist_samples) 63 N.clip(res, 1.0/(nsamples+2), (nsamples+1.0)/(nsamples+2), res) 64 elif self._correction is None: 65 pass 66 else: 67 raise ValueError, \ 68 '%r is incorrect value for correction parameter of %s' \ 69 % (self._correction, self.__class__.__name__) 70 return res

71

72 73 -def _pvalue(x, cdf_func, tail, return_tails=False, name=None):

74 """Helper function to return p-value(x) given cdf and tail 75 76 :Parameters: 77 cdf_func : callable 78 Function to be used to derive cdf values for x 79 tail : str ('left', 'right', 'any', 'both') 80 Which tail of the distribution to report. For 'any' and 'both' 81 it chooses the tail it belongs to based on the comparison to 82 p=0.5. In the case of 'any' significance is taken like in a 83 one-tailed test. 84 return_tails : bool 85 If True, a tuple return (pvalues, tails), where tails contain 86 1s if value was from the right tail, and 0 if the value was 87 from the left tail. 88 """ 89 is_scalar = N.isscalar(x) 90 if is_scalar: 91 x = [x] 92 93 cdf = cdf_func(x) 94 95 if __debug__ and 'CHECK_STABILITY' in debug.active: 96 cdf_min, cdf_max = N.min(cdf), N.max(cdf) 97 if cdf_min < 0 or cdf_max > 1.0: 98 s = ('', ' for %s' % name)[int(name is not None)] 99 warning('Stability check of cdf %s failed%s. Min=%s, max=%s' % \ 100 (cdf_func, s, cdf_min, cdf_max)) 101 102 # no escape but to assure that CDF is in the right range. Some 103 # distributions from scipy tend to jump away from [0,1] 104 cdf = N.clip(cdf, 0, 1.0) 105 106 if tail == 'left': 107 if return_tails: 108 right_tail = N.zeros(cdf.shape, dtype=bool) 109 elif tail == 'right': 110 cdf = 1 - cdf 111 if return_tails: 112 right_tail = N.ones(cdf.shape, dtype=bool) 113 elif tail in ('any', 'both'): 114 right_tail = (cdf >= 0.5) 115 cdf[right_tail] = 1.0 - cdf[right_tail] 116 if tail == 'both': 117 # we need to half the signficance 118 cdf *= 2 119 120 # Assure that NaNs didn't get significant value 121 cdf[N.isnan(x)] = 1.0 122 if is_scalar: res = cdf[0] 123 else: res = cdf 124 125 if return_tails: 126 return (res, right_tail) 127 else: 128 return res

129

130 131 -class NullDist(ClassWithCollections):

132 """Base class for null-hypothesis testing. 133 134 """ 135 136 # Although base class is not benefiting from states, derived 137 # classes do (MCNullDist). For the sake of avoiding multiple 138 # inheritance and associated headache -- let them all be ClassWithCollections, 139 # performance hit should be negligible in most of the scenarios 140 _ATTRIBUTE_COLLECTIONS = ['states'] 141

142 - def __init__(self, tail='both', **kwargs):

143 """Cheap initialization. 144 145 :Parameter: 146 tail: str ('left', 'right', 'any', 'both') 147 Which tail of the distribution to report. For 'any' and 'both' 148 it chooses the tail it belongs to based on the comparison to 149 p=0.5. In the case of 'any' significance is taken like in a 150 one-tailed test. 151 """ 152 ClassWithCollections.__init__(self, **kwargs) 153 154 self._setTail(tail)

155

156 - def __repr__(self, prefixes=[]):

157 return super(NullDist, self).__repr__( 158 prefixes=["tail=%s" % `self.__tail`] + prefixes)

159 160

161 - def _setTail(self, tail):

162 # sanity check 163 if tail not in ('left', 'right', 'any', 'both'): 164 raise ValueError, 'Unknown value "%s" to `tail` argument.' \ 165 % tail 166 self.__tail = tail

167 168

169 - def fit(self, measure, wdata, vdata=None):

170 """Implement to fit the distribution to the data.""" 171 raise NotImplementedError

172 173

174 - def cdf(self, x):

175 """Implementations return the value of the cumulative distribution 176 function (left or right tail dpending on the setting). 177 """ 178 raise NotImplementedError

179 180

181 - def p(self, x, **kwargs):

182 """Returns the p-value for values of `x`. 183 Returned values are determined left, right, or from any tail 184 depending on the constructor setting. 185 186 In case a `FeaturewiseDatasetMeasure` was used to estimate the 187 distribution the method returns an array. In that case `x` can be 188 a scalar value or an array of a matching shape. 189 """ 190 return _pvalue(x, self.cdf, self.__tail, **kwargs)

191 192 193 tail = property(fget=lambda x:x.__tail, fset=_setTail)

194

195 196 -class MCNullDist(NullDist):

197 """Null-hypothesis distribution is estimated from randomly permuted data labels. 198 199 The distribution is estimated by calling fit() with an appropriate 200 `DatasetMeasure` or `TransferError` instance and a training and a 201 validation dataset (in case of a `TransferError`). For a customizable 202 amount of cycles the training data labels are permuted and the 203 corresponding measure computed. In case of a `TransferError` this is the 204 error when predicting the *correct* labels of the validation dataset. 205 206 The distribution can be queried using the `cdf()` method, which can be 207 configured to report probabilities/frequencies from `left` or `right` tail, 208 i.e. fraction of the distribution that is lower or larger than some 209 critical value. 210 211 This class also supports `FeaturewiseDatasetMeasure`. In that case `cdf()` 212 returns an array of featurewise probabilities/frequencies. 213 """ 214 215 _DEV_DOC = """ 216 TODO automagically decide on the number of samples/permutations needed 217 Caution should be paid though since resultant distributions might be 218 quite far from some conventional ones (e.g. Normal) -- it is expected to 219 them to be bimodal (or actually multimodal) in many scenarios. 220 """ 221 222 dist_samples = StateVariable(enabled=False, 223 doc='Samples obtained for each permutation') 224

225 - def __init__(self, dist_class=Nonparametric, permutations=100, **kwargs):

226 """Initialize Monte-Carlo Permutation Null-hypothesis testing 227 228 :Parameters: 229 dist_class: class 230 This can be any class which provides parameters estimate 231 using `fit()` method to initialize the instance, and 232 provides `cdf(x)` method for estimating value of x in CDF. 233 All distributions from SciPy's 'stats' module can be used. 234 permutations: int 235 This many permutations of label will be performed to 236 determine the distribution under the null hypothesis. 237 """ 238 NullDist.__init__(self, **kwargs) 239 240 self._dist_class = dist_class 241 self._dist = [] # actual distributions 242 243 self.__permutations = permutations 244 """Number of permutations to compute the estimate the null 245 distribution."""

246

247 - def __repr__(self, prefixes=[]):

248 prefixes_ = ["permutations=%s" % self.__permutations] 249 if self._dist_class != Nonparametric: 250 prefixes_.insert(0, 'dist_class=%s' % `self._dist_class`) 251 return super(MCNullDist, self).__repr__( 252 prefixes=prefixes_ + prefixes)

253 254

255 - def fit(self, measure, wdata, vdata=None):

256 """Fit the distribution by performing multiple cycles which repeatedly 257 permuted labels in the training dataset. 258 259 :Parameters: 260 measure: (`Featurewise`)`DatasetMeasure` | `TransferError` 261 TransferError instance used to compute all errors. 262 wdata: `Dataset` which gets permuted and used to compute the 263 measure/transfer error multiple times. 264 vdata: `Dataset` used for validation. 265 If provided measure is assumed to be a `TransferError` and 266 working and validation dataset are passed onto it. 267 """ 268 dist_samples = [] 269 """Holds the values for randomized labels.""" 270 271 # Needs to be imported here upon demand due to circular imports 272 # TODO: place MC into a separate module 273 from mvpa.clfs.base import DegenerateInputError, FailedToTrainError 274 275 # decide on the arguments to measure 276 if not vdata is None: 277 measure_args = [vdata, wdata] 278 else: 279 measure_args = [wdata] 280 281 # estimate null-distribution 282 skipped = 0 # # of skipped permutations 283 for p in xrange(self.__permutations): 284 # new permutation all the time 285 # but only permute the training data and keep the testdata constant 286 # 287 if __debug__: 288 debug('STATMC', "Doing %i permutations: %i" \ 289 % (self.__permutations, p+1), cr=True) 290 291 # TODO this really needs to be more clever! If data samples are 292 # shuffled within a class it really makes no difference for the 293 # classifier, hence the number of permutations to estimate the 294 # null-distribution of transfer errors can be reduced dramatically 295 # when the *right* permutations (the ones that matter) are done. 296 wdata.permuteLabels(True, perchunk=False) 297 298 # compute and store the measure of this permutation 299 # assume it has `TransferError` interface 300 try: 301 dist_samples.append(measure(*measure_args)) 302 except (DegenerateInputError, FailedToTrainError), e: 303 if __debug__: 304 debug('STATMC', " skipped", cr=True) 305 warning("Failed to estimate %s on %s, due to %s. " 306 "Permutation %d skipped." % 307 (measure, measure_args, e, p)) 308 skipped += 1 309 continue 310 311 if __debug__: 312 debug('STATMC', ' Skipped: %d permutations' % skipped) 313 314 # restore original labels 315 wdata.permuteLabels(False, perchunk=False) 316 317 # store samples 318 self.dist_samples = dist_samples = N.asarray(dist_samples) 319 320 # fit distribution per each element 321 322 # to decide either it was done on scalars or vectors 323 shape = dist_samples.shape 324 nshape = len(shape) 325 # if just 1 dim, original data was scalar, just create an 326 # artif dimension for it 327 if nshape == 1: 328 dist_samples = dist_samples[:, N.newaxis] 329 330 # fit per each element. 331 # XXX could be more elegant? may be use N.vectorize? 332 dist_samples_rs = dist_samples.reshape((shape[0], -1)) 333 dist = [] 334 for samples in dist_samples_rs.T: 335 params = self._dist_class.fit(samples) 336 if __debug__ and 'STAT' in debug.active: 337 debug('STAT', 'Estimated parameters for the %s are %s' 338 % (self._dist_class, str(params))) 339 dist.append(self._dist_class(*params)) 340 self._dist = dist

341 342

343 - def cdf(self, x):

344 """Return value of the cumulative distribution function at `x`. 345 """ 346 if self._dist is None: 347 # XXX We might not want to descriminate that way since 348 # usually generators also have .cdf where they rely on the 349 # default parameters. But then what about Nonparametric 350 raise RuntimeError, "Distribution has to be fit first" 351 352 is_scalar = N.isscalar(x) 353 if is_scalar: 354 x = [x] 355 356 x = N.asanyarray(x) 357 xshape = x.shape 358 # assure x is a 1D array now 359 x = x.reshape((-1,)) 360 361 if len(self._dist) != len(x): 362 raise ValueError, 'Distribution was fit for structure with %d' \ 363 ' elements, whenever now queried with %d elements' \ 364 % (len(self._dist), len(x)) 365 366 # extract cdf values per each element 367 cdfs = [ dist.cdf(v) for v, dist in zip(x, self._dist) ] 368 return N.array(cdfs).reshape(xshape)

369 370

371 - def clean(self):

372 """Clean stored distributions 373 374 Storing all of the distributions might be too expensive 375 (e.g. in case of Nonparametric), and the scope of the object 376 might be too broad to wait for it to be destroyed. Clean would 377 bind dist_samples to empty list to let gc revoke the memory. 378 """ 379 self._dist = []

380

381 382 383 -class FixedNullDist(NullDist):

384 """Proxy/Adaptor class for SciPy distributions. 385 386 All distributions from SciPy's 'stats' module can be used with this class. 387 388 >>> import numpy as N 389 >>> from scipy import stats 390 >>> from mvpa.clfs.stats import FixedNullDist 391 >>> 392 >>> dist = FixedNullDist(stats.norm(loc=2, scale=4)) 393 >>> dist.p(2) 394 0.5 395 >>> 396 >>> dist.cdf(N.arange(5)) 397 array([ 0.30853754, 0.40129367, 0.5 , 0.59870633, 0.69146246]) 398 >>> 399 >>> dist = FixedNullDist(stats.norm(loc=2, scale=4), tail='right') 400 >>> dist.p(N.arange(5)) 401 array([ 0.69146246, 0.59870633, 0.5 , 0.40129367, 0.30853754]) 402 """

403 - def __init__(self, dist, **kwargs):

404 """ 405 :Parameter: 406 dist: distribution object 407 This can be any object the has a `cdf()` method to report the 408 cumulative distribition function values. 409 """ 410 NullDist.__init__(self, **kwargs) 411 412 self._dist = dist

413 414

415 - def fit(self, measure, wdata, vdata=None):

416 """Does nothing since the distribution is already fixed.""" 417 pass

418 419

420 - def cdf(self, x):

421 """Return value of the cumulative distribution function at `x`. 422 """ 423 return self._dist.cdf(x)

424 425

426 - def __repr__(self, prefixes=[]):

427 prefixes_ = ["dist=%s" % `self._dist`] 428 return super(FixedNullDist, self).__repr__( 429 prefixes=prefixes_ + prefixes)

430

431 432 -class AdaptiveNullDist(FixedNullDist):

433 """Adaptive distribution which adjusts parameters according to the data 434 435 WiP: internal implementation might change 436 """

437 - def fit(self, measure, wdata, vdata=None):

438 """Cares about dimensionality of the feature space in measure 439 """ 440 441 try: 442 nfeatures = len(measure.feature_ids) 443 except ValueError: # XXX 444 nfeatures = N.prod(wdata.shape[1:]) 445 446 dist_gen = self._dist 447 if not hasattr(dist_gen, 'fit'): # frozen already 448 dist_gen = dist_gen.dist # rv_frozen at least has it ;) 449 450 args, kwargs = self._adapt(nfeatures, measure, wdata, vdata) 451 if __debug__: 452 debug('STAT', 'Adapted parameters for %s to be %s, %s' 453 % (dist_gen, args, kwargs)) 454 self._dist = dist_gen(*args, **kwargs)

455 456

457 - def _adapt(self, nfeatures, measure, wdata, vdata=None):

458 raise NotImplementedError

459

460 461 -class AdaptiveRDist(AdaptiveNullDist):

462 """Adaptive rdist: params are (nfeatures-1, 0, 1) 463 """ 464

465 - def _adapt(self, nfeatures, measure, wdata, vdata=None):

466 return (nfeatures-1, 0, 1), {}

467 468 # XXX: RDist has stability issue, just run 469 # python -c "import scipy.stats; print scipy.stats.rdist(541,0,1).cdf(0.72)" 470 # to get some improbable value, so we need to take care about that manually 471 # here

472 - def cdf(self, x):

473 cdf_ = self._dist.cdf(x) 474 bad_values = N.where(N.abs(cdf_)>1) 475 # XXX there might be better implementation (faster/elegant) using N.clip, 476 # the only problem is that instability results might flip the sign 477 # arbitrarily 478 if len(bad_values[0]): 479 # in this distribution we have mean at 0, so we can take that easily 480 # into account 481 cdf_bad = cdf_[bad_values] 482 x_bad = x[bad_values] 483 cdf_bad[x_bad<0] = 0.0 484 cdf_bad[x_bad>=0] = 1.0 485 cdf_[bad_values] = cdf_bad 486 return cdf_

487

488 489 -class AdaptiveNormal(AdaptiveNullDist):

490 """Adaptive Normal Distribution: params are (0, sqrt(1/nfeatures)) 491 """ 492

493 - def _adapt(self, nfeatures, measure, wdata, vdata=None):

494 return (0, 1.0/N.sqrt(nfeatures)), {}

495 496 497 if externals.exists('scipy'): 498 from mvpa.support.stats import scipy 499 from scipy.stats import kstest 500 """ 501 Thoughts: 502 503 So we can use `scipy.stats.kstest` (Kolmogorov-Smirnov test) to 504 check/reject H0 that samples come from a given distribution. But 505 since it is based on a full range of data, we might better of with 506 some ad-hoc checking by the detection power of the values in the 507 tail of a tentative distribution. 508 509 """ 510 511 # We need a way to fixate estimation of some parameters 512 # (e.g. mean) so lets create a simple proxy, or may be class 513 # generator later on, which would take care about punishing change 514 # from the 'right' arguments 515 516 import scipy

517 518 - class rv_semifrozen(object):

519 """Helper proxy-class to fit distribution when some parameters are known 520 521 It is an ugly hack with snippets of code taken from scipy, which is 522 Copyright (c) 2001, 2002 Enthought, Inc. 523 and is distributed under BSD license 524 http://www.scipy.org/License_Compatibility 525 """ 526

527 - def __init__(self, dist, loc=None, scale=None, args=None):

528 529 self._dist = dist 530 # loc and scale 531 theta = (loc, scale) 532 # args 533 Narg_ = dist.numargs 534 if args is not None: 535 Narg = len(args) 536 if Narg > Narg_: 537 raise ValueError, \ 538 'Distribution %s has only %d arguments. Got %d' \ 539 % (dist, Narg_, Narg) 540 args += (None,) * (Narg_ - Narg) 541 else: 542 args = (None,) * Narg_ 543 544 args_i = [i for i,v in enumerate(args) if v is None] 545 self._fargs = (list(args+theta), args_i) 546 """Arguments which should get some fixed value"""

547 548

549 - def __call__(self, *args, **kwargs):

550 """Upon call mimic call to get actual rv_frozen distribution 551 """ 552 return self._dist(*args, **kwargs)

553 554

555 - def nnlf(self, theta, x):

556 # - sum (log pdf(x, theta),axis=0) 557 # where theta are the parameters (including loc and scale) 558 # 559 fargs, fargs_i = self._fargs 560 try: 561 i=-1 562 if fargs[-1] is not None: 563 scale = fargs[-1] 564 else: 565 scale = theta[i] 566 i -= 1 567 568 if fargs[-2] is not None: 569 loc = fargs[-2] 570 else: 571 loc = theta[i] 572 i -= 1 573 574 args = theta[:i+1] 575 # adjust args if there were fixed 576 for i,a in zip(fargs_i, args): 577 fargs[i] = a 578 args = fargs[:-2] 579 580 except IndexError: 581 raise ValueError, "Not enough input arguments." 582 if not self._argcheck(*args) or scale <= 0: 583 return N.inf 584 x = N.asarray((x-loc) / scale) 585 cond0 = (x <= self.a) | (x >= self.b) 586 if (N.any(cond0)): 587 return N.inf 588 else: 589 return self._nnlf(x, *args) + len(x)*N.log(scale)

590

591 - def fit(self, data, *args, **kwds):

592 loc0, scale0 = map(kwds.get, ['loc', 'scale'], [0.0, 1.0]) 593 fargs, fargs_i = self._fargs 594 Narg = len(args) 595 Narg_ = self.numargs 596 if Narg != Narg_: 597 if Narg > Narg_: 598 raise ValueError, "Too many input arguments." 599 else: 600 args += (1.0,)*(self.numargs-Narg) 601 602 # Provide only those args which are not fixed, and 603 # append location and scale (if not fixed) at the end 604 if len(fargs_i) != Narg_: 605 x0 = tuple([args[i] for i in fargs_i]) 606 else: 607 x0 = args 608 if fargs[-2] is None: x0 = x0 + (loc0,) 609 if fargs[-1] is None: x0 = x0 + (scale0,) 610 611 opt_x = scipy.optimize.fmin(self.nnlf, x0, args=(N.ravel(data),), disp=0) 612 613 # reconstruct back 614 i = 0 615 loc, scale = fargs[-2:] 616 if fargs[-1] is None: 617 i -= 1 618 scale = opt_x[i] 619 if fargs[-2] is None: 620 i -= 1 621 loc = opt_x[i] 622 623 # assign those which weren't fixed 624 for i in fargs_i: 625 fargs[i] = opt_x[i] 626 627 #raise ValueError 628 opt_x = N.hstack((fargs[:-2], (loc, scale))) 629 return opt_x

630 631

632 - def __setattr__(self, a, v):

633 if not a in ['_dist', '_fargs', 'fit', 'nnlf']: 634 self._dist.__setattr__(a, v) 635 else: 636 object.__setattr__(self, a, v)

637 638

639 - def __getattribute__(self, a):

640 """We need to redirect all queries correspondingly 641 """ 642 if not a in ['_dist', '_fargs', 'fit', 'nnlf']: 643 return getattr(self._dist, a) 644 else: 645 return object.__getattribute__(self, a)

646

647 648 649 - def matchDistribution(data, nsamples=None, loc=None, scale=None, 650 args=None, test='kstest', distributions=None, 651 **kwargs):

652 """Determine best matching distribution. 653 654 Can be used for 'smelling' the data, as well to choose a 655 parametric distribution for data obtained from non-parametric 656 testing (e.g. `MCNullDist`). 657 658 WiP: use with caution, API might change 659 660 :Parameters: 661 data : N.ndarray 662 Array of the data for which to deduce the distribution. It has 663 to be sufficiently large to make a reliable conclusion 664 nsamples : int or None 665 If None -- use all samples in data to estimate parametric 666 distribution. Otherwise use only specified number randomly selected 667 from data. 668 loc : float or None 669 Loc for the distribution (if known) 670 scale : float or None 671 Scale for the distribution (if known) 672 test : basestring 673 What kind of testing to do. Choices: 674 'p-roc' : detection power for a given ROC. Needs two 675 parameters: `p=0.05` and `tail='both'` 676 'kstest' : 'full-body' distribution comparison. The best 677 choice is made by minimal reported distance after estimating 678 parameters of the distribution. Parameter `p=0.05` sets 679 threshold to reject null-hypothesis that distribution is the 680 same. 681 WARNING: older versions (e.g. 0.5.2 in etch) of scipy have 682 incorrect kstest implementation and do not function 683 properly 684 distributions : None or list of basestring or tuple(basestring, dict) 685 Distributions to check. If None, all known in scipy.stats 686 are tested. If distribution is specified as a tuple, then 687 it must contain name and additional parameters (name, loc, 688 scale, args) in the dictionary. Entry 'scipy' adds all known 689 in scipy.stats. 690 **kwargs 691 Additional arguments which are needed for each particular test 692 (see above) 693 694 :Example: 695 data = N.random.normal(size=(1000,1)); 696 matches = matchDistribution( 697 data, 698 distributions=['rdist', 699 ('rdist', {'name':'rdist_fixed', 700 'loc': 0.0, 701 'args': (10,)})], 702 nsamples=30, test='p-roc', p=0.05) 703 """ 704 705 # Handle parameters 706 _KNOWN_TESTS = ['p-roc', 'kstest'] 707 if not test in _KNOWN_TESTS: 708 raise ValueError, 'Unknown kind of test %s. Known are %s' \ 709 % (test, _KNOWN_TESTS) 710 711 data = N.ravel(data) 712 # data sampled 713 if nsamples is not None: 714 if __debug__: 715 debug('STAT', 'Sampling %d samples from data for the ' \ 716 'estimation of the distributions parameters' % nsamples) 717 indexes_selected = (N.random.sample(nsamples)*len(data)).astype(int) 718 data_selected = data[indexes_selected] 719 else: 720 indexes_selected = N.arange(len(data)) 721 data_selected = data 722 723 p_thr = kwargs.get('p', 0.05) 724 if test == 'p-roc': 725 tail = kwargs.get('tail', 'both') 726 data_p = _pvalue(data, Nonparametric(data).cdf, tail) 727 data_p_thr = N.abs(data_p) <= p_thr 728 true_positives = N.sum(data_p_thr) 729 if true_positives == 0: 730 raise ValueError, "Provided data has no elements in non-" \ 731 "parametric distribution under p<=%.3f. Please " \ 732 "increase the size of data or value of p" % p 733 if __debug__: 734 debug('STAT_', 'Number of positives in non-parametric ' 735 'distribution is %d' % true_positives) 736 737 if distributions is None: 738 distributions = ['scipy'] 739 740 # lets see if 'scipy' entry was in there 741 try: 742 scipy_ind = distributions.index('scipy') 743 distributions.pop(scipy_ind) 744 sp_dists = scipy.stats.distributions.__all__ 745 sp_version = externals.versions['scipy'] 746 if sp_version >= '0.9.0': 747 for d_ in ['ncf']: 748 if d_ in sp_dists: 749 warning("Not considering %s distribution because of " 750 "known issues in scipy %s" % (d_, sp_version)) 751 _ = sp_dists.pop(sp_dists.index(d_)) 752 distributions += sp_dists 753 except ValueError: 754 pass 755 756 results = [] 757 for d in distributions: 758 dist_gen, loc_, scale_, args_ = None, loc, scale, args 759 if isinstance(d, basestring): 760 dist_gen = d 761 dist_name = d 762 elif isinstance(d, tuple): 763 if not (len(d)==2 and isinstance(d[1], dict)): 764 raise ValueError,\ 765 "Tuple specification of distribution must be " \ 766 "(d, {params}). Got %s" % (d,) 767 dist_gen = d[0] 768 loc_ = d[1].get('loc', loc) 769 scale_ = d[1].get('scale', scale) 770 args_ = d[1].get('args', args) 771 dist_name = d[1].get('name', str(dist_gen)) 772 else: 773 dist_gen = d 774 dist_name = str(d) 775 776 # perform actions which might puke for some distributions 777 try: 778 dist_gen_ = getattr(scipy.stats, dist_gen) 779 # specify distribution 'optimizer'. If loc or scale was provided, 780 # use home-brewed rv_semifrozen 781 if args_ is not None or loc_ is not None or scale_ is not None: 782 dist_opt = rv_semifrozen(dist_gen_, loc=loc_, scale=scale_, args=args_) 783 else: 784 dist_opt = dist_gen_ 785 dist_params = dist_opt.fit(data_selected) 786 if __debug__: 787 debug('STAT__', 788 'Got distribution parameters %s for %s' 789 % (dist_params, dist_name)) 790 if test == 'p-roc': 791 cdf_func = lambda x: dist_gen_.cdf(x, *dist_params) 792 # We need to compare detection under given p 793 cdf_p = N.abs(_pvalue(data, cdf_func, tail, name=dist_gen)) 794 cdf_p_thr = cdf_p <= p_thr 795 D, p = N.sum(N.abs(data_p_thr - cdf_p_thr))*1.0/true_positives, 1 796 if __debug__: res_sum = 'D=%.2f' % D 797 elif test == 'kstest': 798 D, p = kstest(data, d, args=dist_params) 799 if __debug__: res_sum = 'D=%.3f p=%.3f' % (D, p) 800 except (TypeError, ValueError, AttributeError, 801 NotImplementedError), e:#Exception, e: 802 if __debug__: 803 debug('STAT__', 804 'Testing for %s distribution failed due to %s' 805 % (d, str(e))) 806 continue 807 808 if p > p_thr and not N.isnan(D): 809 results += [ (D, dist_gen, dist_name, dist_params) ] 810 if __debug__: 811 debug('STAT_', 'Testing for %s dist.: %s' % (dist_name, res_sum)) 812 else: 813 if __debug__: 814 debug('STAT__', 'Cannot consider %s dist. with %s' 815 % (d, res_sum)) 816 continue 817 818 # sort in ascending order, so smaller is better 819 results.sort() 820 821 if __debug__ and 'STAT' in debug.active: 822 # find the best and report it 823 nresults = len(results) 824 sresult = lambda r:'%s(%s)=%.2f' % (r[1], ', '.join(map(str, r[3])), r[0]) 825 if nresults>0: 826 nnextbest = min(2, nresults-1) 827 snextbest = ', '.join(map(sresult, results[1:1+nnextbest])) 828 debug('STAT', 'Best distribution %s. Next best: %s' 829 % (sresult(results[0]), snextbest)) 830 else: 831 debug('STAT', 'Could not find suitable distribution') 832 833 # return all the results 834 return results

835 836 837 if externals.exists('pylab'): 838 import pylab as P

839 840 - def plotDistributionMatches(data, matches, nbins=31, nbest=5, 841 expand_tails=8, legend=2, plot_cdf=True, 842 p=None, tail='both'):

843 """Plot best matching distributions 844 845 :Parameters: 846 data : N.ndarray 847 Data which was used to obtain the matches 848 matches : list of tuples 849 Sorted matches as provided by matchDistribution 850 nbins : int 851 Number of bins in the histogram 852 nbest : int 853 Number of top matches to plot 854 expand_tails : int 855 How many bins away to add to parametrized distributions 856 plots 857 legend : int 858 Either to provide legend and statistics in the legend. 859 1 -- just lists distributions. 860 2 -- adds distance measure 861 3 -- tp/fp/fn in the case if p is provided 862 plot_cdf : bool 863 Either to plot cdf for data using non-parametric distribution 864 p : float or None 865 If not None, visualize null-hypothesis testing (given p). 866 Bars in the histogram which fall under given p are colored 867 in red. False positives and false negatives are marked as 868 triangle up and down symbols correspondingly 869 tail : ('left', 'right', 'any', 'both') 870 If p is not None, the choise of tail for null-hypothesis 871 testing 872 873 :Returns: tuple(histogram, list of lines) 874 """ 875 876 hist = P.hist(data, nbins, normed=1, align='center') 877 data_range = [N.min(data), N.max(data)] 878 879 # x's 880 x = hist[1] 881 dx = x[expand_tails] - x[0] # how much to expand tails by 882 x = N.hstack((x[:expand_tails] - dx, x, x[-expand_tails:] + dx)) 883 884 nonparam = Nonparametric(data) 885 # plot cdf 886 if plot_cdf: 887 P.plot(x, nonparam.cdf(x), 'k--', linewidth=1) 888 889 p_thr = p 890 891 data_p = _pvalue(data, nonparam.cdf, tail) 892 data_p_thr = (data_p <= p_thr).ravel() 893 894 x_p = _pvalue(x, Nonparametric(data).cdf, tail) 895 x_p_thr = N.abs(x_p) <= p_thr 896 # color bars which pass thresholding in red 897 for thr, bar in zip(x_p_thr[expand_tails:], hist[2]): 898 bar.set_facecolor(('w','r')[int(thr)]) 899 900 if not len(matches): 901 # no matches were provided 902 warning("No matching distributions were provided -- nothing to plot") 903 return (hist, ) 904 905 lines = [] 906 labels = [] 907 for i in xrange(min(nbest, len(matches))): 908 D, dist_gen, dist_name, params = matches[i] 909 dist = getattr(scipy.stats, dist_gen)(*params) 910 911 label = '%s' % (dist_name) 912 if legend > 1: label += '(D=%.2f)' % (D) 913 914 xcdf_p = N.abs(_pvalue(x, dist.cdf, tail)) 915 xcdf_p_thr = (xcdf_p <= p_thr).ravel() 916 917 if p is not None and legend > 2: 918 # We need to compare detection under given p 919 data_cdf_p = N.abs(_pvalue(data, dist.cdf, tail)) 920 data_cdf_p_thr = (data_cdf_p <= p_thr).ravel() 921 922 # true positives 923 tp = N.logical_and(data_cdf_p_thr, data_p_thr) 924 # false positives 925 fp = N.logical_and(data_cdf_p_thr, ~data_p_thr) 926 # false negatives 927 fn = N.logical_and(~data_cdf_p_thr, data_p_thr) 928 929 label += ' tp/fp/fn=%d/%d/%d)' % \ 930 tuple(map(N.sum, [tp,fp,fn])) 931 932 pdf = dist.pdf(x) 933 line = P.plot(x, pdf, '-', linewidth=2, label=label) 934 color = line[0].get_color() 935 936 if plot_cdf: 937 cdf = dist.cdf(x) 938 P.plot(x, cdf, ':', linewidth=1, color=color, label=label) 939 940 # TODO: decide on tp/fp/fn by not centers of the bins but 941 # by the values in data in the ranges covered by 942 # those bins. Then it would correspond to the values 943 # mentioned in the legend 944 if p is not None: 945 # true positives 946 xtp = N.logical_and(xcdf_p_thr, x_p_thr) 947 # false positives 948 xfp = N.logical_and(xcdf_p_thr, ~x_p_thr) 949 # false negatives 950 xfn = N.logical_and(~xcdf_p_thr, x_p_thr) 951 952 # no need to plot tp explicitely -- marked by color of the bar 953 # P.plot(x[xtp], pdf[xtp], 'o', color=color) 954 P.plot(x[xfp], pdf[xfp], '^', color=color) 955 P.plot(x[xfn], pdf[xfn], 'v', color=color) 956 957 lines.append(line) 958 labels.append(label) 959 960 if legend: 961 P.legend(lines, labels) 962 963 return (hist, lines)

964

965 #if True: 966 # data = N.random.normal(size=(1000,1)); 967 # matches = matchDistribution( 968 # data, 969 # distributions=['scipy', 970 # ('norm', {'name':'norm_known', 971 # 'scale': 1.0, 972 # 'loc': 0.0})], 973 # nsamples=30, test='p-roc', p=0.05) 974 # P.figure(); plotDistributionMatches(data, matches, nbins=101, 975 # p=0.05, legend=4, nbest=5) 976 977 978 -def autoNullDist(dist):

979 """Cheater for human beings -- wraps `dist` if needed with some 980 NullDist 981 982 tail and other arguments are assumed to be default as in 983 NullDist/MCNullDist 984 """ 985 if dist is None or isinstance(dist, NullDist): 986 return dist 987 elif hasattr(dist, 'fit'): 988 if __debug__: 989 debug('STAT', 'Wrapping %s into MCNullDist' % dist) 990 return MCNullDist(dist) 991 else: 992 if __debug__: 993 debug('STAT', 'Wrapping %s into FixedNullDist' % dist) 994 return FixedNullDist(dist)

995

996 997 # if no scipy, we need nanmean 998 -def _chk_asarray(a, axis):

999 if axis is None: 1000 a = N.ravel(a) 1001 outaxis = 0 1002 else: 1003 a = N.asarray(a) 1004 outaxis = axis 1005 return a, outaxis

1006

1007 -def nanmean(x, axis=0):

1008 """Compute the mean over the given axis ignoring nans. 1009 1010 :Parameters: 1011 x : ndarray 1012 input array 1013 axis : int 1014 axis along which the mean is computed. 1015 1016 :Results: 1017 m : float 1018 the mean.""" 1019 x, axis = _chk_asarray(x,axis) 1020 x = x.copy() 1021 Norig = x.shape[axis] 1022 factor = 1.0-N.sum(N.isnan(x),axis)*1.0/Norig 1023 1024 x[N.isnan(x)] = 0 1025 return N.mean(x,axis)/factor

1026

Source Code for Module mvpa.clfs.stats