Package mvpa :: Package support :: Module stats

[hide private]
[frames] | no frames]

# 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