# [Scipy-svn] r5167 - in trunk/scipy/stats: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Nov 22 10:10:19 CST 2008

```Author: josef
Date: 2008-11-22 10:10:16 -0600 (Sat, 22 Nov 2008)
New Revision: 5167

Modified:
trunk/scipy/stats/distributions.py
trunk/scipy/stats/tests/test_continuous_basic.py
Log:
genextreme: correct _argcheck for c close to zero, allow array argument in _cdf

Modified: trunk/scipy/stats/distributions.py
===================================================================
--- trunk/scipy/stats/distributions.py	2008-11-22 14:40:43 UTC (rev 5166)
+++ trunk/scipy/stats/distributions.py	2008-11-22 16:10:16 UTC (rev 5167)
@@ -14,7 +14,7 @@
from numpy import alltrue, where, arange, put, putmask, \
ravel, take, ones, sum, shape, product, repeat, reshape, \
zeros, floor, logical_and, log, sqrt, exp, arctanh, tan, sin, arcsin, \
-     arctan, tanh, ndarray, cos, cosh, sinh, newaxis, array
+     arctan, tanh, ndarray, cos, cosh, sinh, newaxis, array, log1p, expm1
from numpy import atleast_1d, polyval, angle, ceil, place, extract, \
any, argsort, argmax, vectorize, r_, asarray, nan, inf, pi, isnan, isinf, \
power
@@ -49,6 +49,8 @@
'poisson', 'planck', 'boltzmann', 'randint', 'zipf', 'dlaplace',
]

+floatinfo = numpy.finfo(float)
+
errp = special.errprint
arr = asarray
gam = special.gamma
@@ -1744,7 +1746,7 @@

## Generalized Extreme Value
##  c=0 is just gumbel distribution.
-##  This version does not accept c==0
+##  This version does now accept c==0
##  Use gumbel_r for c==0

# new version by Per Brodtkorb, see ticket:767
@@ -1753,18 +1755,23 @@

class genextreme_gen(rv_continuous):
def _argcheck(self, c):
-        self.b = where(c > 0, 1.0 / c, inf)
-        self.a = where(c < 0, 1.0 / c, -inf)
-        return (c==c) #True #(c!=0) #see ticket:793
+        min = np.minimum
+        max = np.maximum
+        sml = floatinfo.machar.xmin
+        #self.b = where(c > 0, 1.0 / c,inf)
+        #self.a = where(c < 0, 1.0 / c, -inf)
+        self.b = where(c > 0, 1.0 / max(c, sml),inf)
+        self.a = where(c < 0, 1.0 / min(c,-sml), -inf)
+        return where(abs(c)==inf, 0, 1) #True #(c!=0)
def _pdf(self, x, c):
##        ex2 = 1-c*x
##        pex2 = pow(ex2,1.0/c)
##        p2 = exp(-pex2)*pex2/ex2
##        return p2
cx = c*x
-        # Note: fit method requires that _pdf accepts vector x
-        logex2 = where((x==x)*(c==0),0.0,np.log1p(-cx))
-        logpex2 = where((x==x)*(c==0),-x,logex2/c)
+
+        logex2 = where((c==0)*(x==x),0.0,log1p(-cx))
+        logpex2 = where((c==0)*(x==x),-x,logex2/c)
pex2 = exp(logpex2)
# % Handle special cases
logpdf = where((cx==1) | (cx==-inf),-inf,-pex2+logpex2-logex2)
@@ -1775,15 +1782,13 @@

def _cdf(self, x, c):
#return exp(-pow(1-c*x,1.0/c))
-        loglogcdf = where(c==0,-x,np.log1p(-c*x)/c)
+        loglogcdf = where((c==0)*(x==x),-x,log1p(-c*x)/c)
return exp(-exp(loglogcdf))

def _ppf(self, q, c):
#return 1.0/c*(1.-(-log(q))**c)
-        logq = log(q);
-        x = -log(-logq)
-        # _rvs requires that _ppf allows vectorized q
-        return where((q==q)*(c==0),x,-np.expm1(-c*x)/c)
+        x = -log(-log(q))
+        return where((c==0)*(x==x),x,-expm1(-c*x)/c)
def _stats(self,c):

g = lambda n : gam(n*c+1)
@@ -1792,9 +1797,9 @@
g3 = g(3);
g4 = g(4)
g2mg12 = where(abs(c)<1e-7,(c*pi)**2.0/6.0,g2-g1**2.0)
-        gam2k = where(abs(c)<1e-7,pi**2.0/6.0, np.expm1(gamln(2.0*c+1.0)-2*gamln(c+1.0))/c**2.0);
+        gam2k = where(abs(c)<1e-7,pi**2.0/6.0, expm1(gamln(2.0*c+1.0)-2*gamln(c+1.0))/c**2.0);
eps = 1e-14
-        gamk = where(abs(c)<eps,-_EULER,np.expm1(gamln(c+1))/c)
+        gamk = where(abs(c)<eps,-_EULER,expm1(gamln(c+1))/c)

m = where(c<-1.0,nan,-gamk)
v = where(c<-0.5,nan,g1**2.0*gam2k)

Modified: trunk/scipy/stats/tests/test_continuous_basic.py
===================================================================
--- trunk/scipy/stats/tests/test_continuous_basic.py	2008-11-22 14:40:43 UTC (rev 5166)
+++ trunk/scipy/stats/tests/test_continuous_basic.py	2008-11-22 16:10:16 UTC (rev 5167)
@@ -117,7 +117,9 @@
##distcont = [
##    ['genextreme', (3.3184017469423535,)],
##    ['genextreme', (0.01,)],
-##    ['genextreme', (0.00001,)]
+##    ['genextreme', (0.00001,)],
+##    ['genextreme', (0.0,)],
+##    ['genextreme', (-0.01,)]
##    ]

def test_cont_basic():

```