Source Code for Module mvpa.support.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  """Fixer for rdist in scipy 
 10  """ 
 11   
 12  __docformat__ = 'restructuredtext' 
 13   
 14  from mvpa.base import externals, warning, cfg 
 15   
 16  if __debug__: 
 17      from mvpa.base import debug 
 18   
 19  if externals.exists('scipy', raiseException=True): 
 20      import scipy 
 21      import scipy.stats 
 22      import scipy.stats as stats 
 23   
 24  if not externals.exists('good scipy.stats.rdist'): 
 25      if __debug__: 
 26          debug("EXT", "Fixing up scipy.stats.rdist") 
 27      # Lets fix it up, future imports of scipy.stats should carry fixed 
 28      # version, isn't python is \emph{evil} ;-) 
 29      import numpy as N 
 30   
 31      from scipy.stats.distributions import rv_continuous 
 32      from scipy import special 
 33      import scipy.integrate 
 34   
 35      # NB: Following function is copied from scipy SVN rev.5236 
 36      #     and fixed with pow -> N.power (thanks Josef!) 
 37      # FIXME: PPF does not work. 
 38 -    class rdist_gen(rv_continuous): 
 39 -        def _pdf(self, x, c): 
 40              return N.power((1.0-x*x),c/2.0-1) / special.beta(0.5,c/2.0) 
 41 -        def _cdf_skip(self, x, c): 
 42              #error inspecial.hyp2f1 for some values see tickets 758, 759 
 43              return 0.5 + x/special.beta(0.5,c/2.0)* \ 
 44                     special.hyp2f1(0.5,1.0-c/2.0,1.5,x*x) 
 45 -        def _munp(self, n, c): 
 46              return (1-(n % 2))*special.beta((n+1.0)/2,c/2.0) 
 47   
 48      # Lets try to avoid at least some of the numerical problems by removing points 
 49      # around edges 
 50      rdist = rdist_gen(a=-1.0, b=1.0, name="rdist", longname="An R-distributed", 
 51                        shapes="c", extradoc=""" 
 52   
 53      R-distribution 
 54   
 55      rdist.pdf(x,c) = (1-x**2)**(c/2-1) / B(1/2, c/2) 
 56      for -1 <= x <= 1, c > 0. 
 57      """ 
 58                        ) 
 59      # Fix up number of arguments for veccdf's vectorize 
 60      if rdist.veccdf.nin == 1: 
 61          if __debug__: 
 62              debug("EXT", "Fixing up veccdf.nin to make 2 for rdist") 
 63          rdist.veccdf.nin = 2 
 64   
 65      scipy.stats.distributions.rdist_gen = scipy.stats.rdist_gen = rdist_gen 
 66      scipy.stats.distributions.rdist = scipy.stats.rdist = rdist 
 67   
 68      try: # Retest 
 69          externals.exists('good scipy.stats.rdist', force=True, 
 70                           raiseException=True) 
 71      except RuntimeError: 
 72          warning("scipy.stats.rdist was not fixed with a monkey-patch. " 
 73                  "It remains broken") 
 74      # Revert so if configuration stored, we know the true flow of things ;) 
 75      cfg.set('externals', 'have good scipy.stats.rdist', 'no') 
 76   
 77   
 78  if not externals.exists('good scipy.stats.rv_discrete.ppf'): 
 79      # Local rebindings for ppf7 (7 is for the scipy version from 
 80      # which code was borrowed) 
 81      arr = N.asarray 
 82      from scipy.stats.distributions import valarray, argsreduce 
 83      from numpy import shape, place, any 
 84   
 85 -    def ppf7(self,q,*args,**kwds): 
 86          """ 
 87          Percent point function (inverse of cdf) at q of the given RV 
 88   
 89          :Parameters: 
 90            q : array-like 
 91                lower tail probability 
 92            arg1, arg2, arg3,... : array-like 
 93                The shape parameter(s) for the distribution (see docstring of the 
 94                instance object for more information) 
 95            loc : array-like, optional 
 96                location parameter (default=0) 
 97   
 98          :Returns: 
 99            k : array-like 
100                quantile corresponding to the lower tail probability, q. 
101   
102          """ 
103          loc = kwds.get('loc') 
104          args, loc = self._rv_discrete__fix_loc(args, loc) 
105          q,loc  = map(arr,(q,loc)) 
106          args = tuple(map(arr,args)) 
107          cond0 = self._argcheck(*args) & (loc == loc) 
108          cond1 = (q > 0) & (q < 1) 
109          cond2 = (q==1) & cond0 
110          cond = cond0 & cond1 
111          output = valarray(shape(cond),value=self.badvalue,typecode='d') 
112          #output type 'd' to handle nin and inf 
113          place(output,(q==0)*(cond==cond), self.a-1) 
114          place(output,cond2,self.b) 
115          if any(cond): 
116              goodargs = argsreduce(cond, *((q,)+args+(loc,))) 
117              loc, goodargs = goodargs[-1], goodargs[:-1] 
118              place(output,cond,self._ppf(*goodargs) + loc) 
119   
120          if output.ndim == 0: 
121              return output[()] 
122          return output 
123   
124      scipy.stats.distributions.rv_discrete.ppf = ppf7 
125      try: 
126          externals.exists('good scipy.stats.rv_discrete.ppf', force=True, 
127                           raiseException=True) 
128      except RuntimeError: 
129          warning("rv_discrete.ppf was not fixed with a monkey-patch. " 
130                  "It remains broken") 
131      cfg.set('externals', 'have good scipy.stats.rv_discrete.ppf', 'no') 
132