Source Code for Module mvpa.support.stats

 1  #emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- 
 2  #ex: set 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 
15  externals.exists('scipy', raiseException=True) 
16   
17  if __debug__: 
18      from mvpa.base import debug 
19   
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 probably due to already mentioned FIXME it is buggy, since if x is close to self.a, 
37      #     squaring it makes it self.a^2, and actually just 1.0 in rdist, so 1-x*x becomes 0 
38      #     which is not raisable to negative non-round powers... 
39      #     as a fix I am initializing .a/.b with values which should avoid such situation 
40      # FIXME: PPF does not work. 
41 -    class rdist_gen(rv_continuous): 
42 -        def _pdf(self, x, c): 
43              return pow((1.0-x*x),c/2.0-1) / special.beta(0.5,c/2.0) 
44 -        def _cdf_skip(self, x, c): 
45              #error inspecial.hyp2f1 for some values see tickets 758, 759 
46              return 0.5 + x/special.beta(0.5,c/2.0)* \ 
47                     special.hyp2f1(0.5,1.0-c/2.0,1.5,x*x) 
48 -        def _munp(self, n, c): 
49              return (1-(n % 2))*special.beta((n+1.0)/2,c/2.0) 
50   
51      # Lets try to avoid at least some of the numerical problems by removing points 
52      # around edges 
53      __eps = N.sqrt(N.finfo(float).eps) 
54      rdist = rdist_gen(a=-1.0+__eps, b=1.0-__eps, name="rdist", longname="An R-distributed", 
55                        shapes="c", extradoc=""" 
56   
57      R-distribution 
58   
59      rdist.pdf(x,c) = (1-x**2)**(c/2-1) / B(1/2, c/2) 
60      for -1 <= x <= 1, c > 0. 
61      """ 
62                        ) 
63      # Fix up number of arguments for veccdf's vectorize 
64      if rdist.veccdf.nin == 1: 
65          if __debug__: 
66              debug("EXT", "Fixing up veccdf.nin to make 2 for rdist") 
67          rdist.veccdf.nin = 2 
68   
69      scipy.stats.distributions.rdist_gen = scipy.stats.rdist_gen = rdist_gen 
70      scipy.stats.distributions.rdist = scipy.stats.rdist = rdist 
71