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

scipy-svn@scip... scipy-svn@scip...
Thu Feb 26 14:24:12 CST 2009

```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",
@@ -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",

@@ -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',
@@ -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

```