Source Code for Module mvpa.measures.anova

  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  """FeaturewiseDatasetMeasure performing a univariate ANOVA.""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13  import numpy as N 
 14   
 15  from mvpa.measures.base import FeaturewiseDatasetMeasure 
 16   
 17  # TODO: Extend with access to functionality from scipy.stats? 
 18  # For binary: 
 19  #  2-sample kolmogorov-smirnof might be interesting 
 20  #   (scipy.stats.ks_2samp) to judge if two conditions are derived 
 21  #   from different distributions (take it as 'activity' vs 'rest'), 
 22  # 
 23  # For binary+multiclass: 
 24  #  kruskal-wallis H-test (scipy.stats.kruskal) 
 25  # 
 26  # and may be some others 
 27   
 28 -class OneWayAnova(FeaturewiseDatasetMeasure): 
 29      """`FeaturewiseDatasetMeasure` that performs a univariate ANOVA. 
 30   
 31      F-scores are computed for each feature as the standard fraction of between 
 32      and within group variances. Groups are defined by samples with unique 
 33      labels. 
 34   
 35      No statistical testing is performed, but raw F-scores are returned as a 
 36      sensitivity map. As usual F-scores have a range of [0,inf] with greater 
 37      values indicating higher sensitivity. 
 38      """ 
 39   
 40 -    def _call(self, dataset, labels=None): 
 41          # This code is based on SciPy's stats.f_oneway() 
 42          # Copyright (c) Gary Strangman.  All rights reserved 
 43          # License: BSD 
 44          # 
 45          # However, it got tweaked and optimized to better fit into PyMVPA. 
 46   
 47          # number of groups 
 48          if labels is None: 
 49              labels = dataset.labels 
 50   
 51          ul = N.unique(labels) 
 52   
 53          na = len(ul) 
 54          bign = float(dataset.nsamples) 
 55          alldata = dataset.samples 
 56   
 57          # total squares of sums 
 58          sostot = N.sum(alldata, axis=0) 
 59          sostot *= sostot 
 60          sostot /= bign 
 61   
 62          # total sum of squares 
 63          sstot = N.sum(alldata * alldata, axis=0) - sostot 
 64   
 65          # between group sum of squares 
 66          ssbn = 0 
 67          for l in ul: 
 68              # all samples for the respective label 
 69              d = alldata[labels == l] 
 70              sos = N.sum(d, axis=0) 
 71              sos *= sos 
 72              ssbn += sos / float(len(d)) 
 73   
 74          ssbn -= sostot 
 75          # within 
 76          sswn = sstot - ssbn 
 77   
 78          # degrees of freedom 
 79          dfbn = na-1 
 80          dfwn = bign - na 
 81   
 82          # mean sums of squares 
 83          msb = ssbn / float(dfbn) 
 84          msw = sswn / float(dfwn) 
 85          f = msb / msw 
 86          # assure no NaNs -- otherwise it leads instead of 
 87          # sane unittest failure (check of NaNs) to crazy 
 88          #   File "mtrand.pyx", line 1661, in mtrand.shuffle 
 89          #  TypeError: object of type 'numpy.int64' has no len() 
 90          # without any sane backtrace 
 91          f[N.isnan(f)] = 0 
 92   
 93          return f 
 94   
 95          # XXX maybe also compute p-values? 
 96          #prob = scipy.stats.fprob(dfbn, dfwn, f) 
 97          #return prob 
 98   
 99   
100 -class CompoundOneWayAnova(OneWayAnova): 
101      """Compound comparisons via univariate ANOVA. 
102   
103      Provides F-scores per each label if compared to the other labels. 
104      """ 
105   
106 -    def _call(self, dataset): 
107          """Computes featurewise f-scores using compound comparisons.""" 
108   
109          orig_labels = dataset.labels 
110          labels = orig_labels.copy() 
111   
112          results = [] 
113          for ul in dataset.uniquelabels: 
114              labels[orig_labels == ul] = 1 
115              labels[orig_labels != ul] = 2 
116              results.append(OneWayAnova._call(self, dataset, labels)) 
117   
118          # features x labels 
119          return N.array(results).T 
120