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

[scipy-svn at scipy.org]

Author: josef
Date: 2009-02-26 13:04:38 -0600 (Thu, 26 Feb 2009)
New Revision: 5599

Modified:
   trunk/scipy/stats/distributions.py
   trunk/scipy/stats/tests/test_distributions.py
Log:
improve stats.expon and exponpow, #807, thanks to Per Brodtkorb and David Huard

Modified: trunk/scipy/stats/distributions.py
===================================================================
--- trunk/scipy/stats/distributions.py	2009-02-26 19:01:09 UTC (rev 5598)
+++ trunk/scipy/stats/distributions.py	2009-02-26 19:04:38 UTC (rev 5599)
@@ -1526,7 +1526,11 @@
     def _pdf(self, x):
         return exp(-x)
     def _cdf(self, x):
-        return 1.0-exp(-x)
+        return -expm1(-x)
+    def _sf(self,x):
+        return exp(-x)
+    def _isf(self,q):
+        return -log(q)
     def _ppf(self, q):
         return -log(1.0-q)
     def _stats(self):
@@ -1576,7 +1580,12 @@
         return exp(1)*b*xbm1 * exp(xb - exp(xb))
     def _cdf(self, x, b):
         xb = arr(x**b)
-        return 1.0-exp(1-exp(xb))
+        return -expm1(-expm1(xb))
+    def _sf(self, x, b):
+        xb = arr(x**b)
+        return exp(-expm1(xb))
+    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)
 exponpow = exponpow_gen(a=0.0,name='exponpow',longname="An exponential power",

Modified: trunk/scipy/stats/tests/test_distributions.py
===================================================================
--- trunk/scipy/stats/tests/test_distributions.py	2009-02-26 19:01:09 UTC (rev 5598)
+++ trunk/scipy/stats/tests/test_distributions.py	2009-02-26 19:04:38 UTC (rev 5599)
@@ -238,6 +238,10 @@
 class TestExpon(TestCase):
     def test_zero(self):
         assert_equal(stats.expon.pdf(0),1)
+        
+    def test_tail(self):  # Regression test for ticket 807
+        assert_equal(stats.expon.cdf(1e-18),  1e-18)
+        assert_equal(stats.expon.isf(stats.expon.sf(40)),  40)
 
 class TestGenExpon(TestCase):
     def test_pdf_unity_area(self):
@@ -251,6 +255,11 @@
         # CDF should always be positive
         cdf = stats.genexpon.cdf(numpy.arange(0, 10, 0.01), 0.5, 0.5, 2.0)
         assert(numpy.all((0 <= cdf) & (cdf <= 1)))
+        
+class TestExponpow(TestCase):
+    def test_tail(self):
+        assert_almost_equal(stats.exponpow.cdf(1e-10,  2.),  1e-20)
+        assert_almost_equal(stats.exponpow.isf(stats.exponpow.sf(5, .8), .8),  5)
 
 class TestDocstring(TestCase):
     def test_docstrings(self):