[Scipy-svn] r5600 - trunk/scipy/stats

[scipy-svn at scipy.org]

Author: josef
Date: 2009-02-26 14:24:09 -0600 (Thu, 26 Feb 2009)
New Revision: 5600

Modified:
   trunk/scipy/stats/distributions.py
Log:
improve numerical precision of distributions using log1p, ex1p, thanks to Per Brodtkorb

Modified: trunk/scipy/stats/distributions.py
===================================================================
--- trunk/scipy/stats/distributions.py	2009-02-26 19:04:38 UTC (rev 5599)
+++ trunk/scipy/stats/distributions.py	2009-02-26 20:24:09 UTC (rev 5600)
@@ -1527,6 +1527,8 @@
         return exp(-x)
     def _cdf(self, x):
         return -expm1(-x)
+    def _ppf(self, q):
+        return -log1p(-q)
     def _sf(self,x):
         return exp(-x)
     def _isf(self,q):
@@ -1556,10 +1558,10 @@
         exc = exp(-x**c)
         return a*c*(1-exc)**arr(a-1) * exc * x**arr(c-1)
     def _cdf(self, x, a, c):
-        exc = exp(-x**c)
-        return arr((1-exc)**a)
+        exm1c = -expm1(-x**c)
+        return arr((exm1c)**a)
     def _ppf(self, q, a, c):
-        return (-log(1-q**(1.0/a)))**arr(1.0/c)
+        return (-log1p(-q**(1.0/a)))**arr(1.0/c)
 exponweib = exponweib_gen(a=0.0,name='exponweib',
                           longname="An exponentiated Weibull",
                           shapes="a,c",extradoc="""
@@ -1587,7 +1589,7 @@
     def _isf(self, x, b):
         return (log1p(-log(x)))**(1./b)
     def _ppf(self, q, b):
-        return pow(log(1.0-log(1.0-q)), 1.0/b)
+        return pow(log1p(-log1p(-q)), 1.0/b)
 exponpow = exponpow_gen(a=0.0,name='exponpow',longname="An exponential power",
                         shapes='b',extradoc="""
 
@@ -1742,9 +1744,9 @@
     def _pdf(self, x, c):
         return c*pow(x,c-1)*exp(-pow(x,c))
     def _cdf(self, x, c):
-        return 1-exp(-pow(x,c))
+        return -expm1(-pow(x,c))
     def _ppf(self, q, c):
-        return pow(-log(1-q),1.0/c)
+        return pow(-log1p(-q),1.0/c)
     def _munp(self, n, c):
         return special.gamma(1.0+n*1.0/c)
     def _entropy(self, c):
@@ -1875,9 +1877,9 @@
 
 class genexpon_gen(rv_continuous):
     def _pdf(self, x, a, b, c):
-        return (a+b*(1-exp(-c*x)))*exp((-a-b)*x+b*(1-exp(-c*x))/c)
+        return (a+b*(-expm1(-c*x)))*exp((-a-b)*x+b*(-expm1(-c*x))/c)
     def _cdf(self, x, a, b, c):
-        return 1.0-exp((-a-b)*x + b*(1-exp(-c*x))/c)
+        return -expm1((-a-b)*x + b*(-expm1(-c*x))/c)
 genexpon = genexpon_gen(a=0.0,name='genexpon',
                         longname='A generalized exponential',
                         shapes='a,b,c',extradoc="""
@@ -3247,9 +3249,9 @@
         #wrong answer with formula, same as in continuous.pdf
         #return gam(n+1)-special.gammainc(1+n,b)
         if n == 1:
-            return (1-(b+1)*np.exp(-b))/(1-np.exp(-b))
+            return (1-(b+1)*exp(-b))/(-expm1(-b))
         elif n == 2:
-            return 2*(1-0.5*(b*b+2*b+2)*np.exp(-b))/(1-np.exp(-b))
+            return 2*(1-0.5*(b*b+2*b+2)*exp(-b))/(-expm1(-b))
         else:
             #return generic for higher moments
             #return rv_continuous._mom1_sc(self,n, b)
@@ -4551,11 +4553,11 @@
         k = floor(x)
         return 1-exp(-lambda_*(k+1))
     def _ppf(self, q, lambda_):
-        val = ceil(-1.0/lambda_ * log(1-q)-1)
+        val = ceil(-1.0/lambda_ * log1p(-q)-1)
         return val
     def _stats(self, lambda_):
         mu = 1/(exp(lambda_)-1)
-        var = exp(-lambda_)/(1-exp(-lambda_))**2
+        var = exp(-lambda_)/(expm1(-lambda_))**2
         g1 = 2*cosh(lambda_/2.0)
         g2 = 4+2*cosh(lambda_)
         return mu, var, g1, g2