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

[scipy-svn at scipy.org]

Author: stefan
Date: 2008-11-03 14:38:44 -0600 (Mon, 03 Nov 2008)
New Revision: 4972

Modified:
   trunk/scipy/stats/distributions.py
   trunk/scipy/stats/tests/test_distributions.py
Log:
Fix generalised bivariate exponential distribution [patch by Konrad Blum].

Modified: trunk/scipy/stats/distributions.py
===================================================================
--- trunk/scipy/stats/distributions.py	2008-11-03 19:51:35 UTC (rev 4971)
+++ trunk/scipy/stats/distributions.py	2008-11-03 20:38:44 UTC (rev 4972)
@@ -1678,17 +1678,24 @@
 
 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*(1-exp(-c*x)))*exp((-a-b)*x+b*(1-exp(-c*x))/c)
     def _cdf(self, x, a, b, c):
-        return 1.0-exp((a-b)*x + b*(1-exp(-c*x))/c)
+        return 1.0-exp((-a-b)*x + b*(1-exp(-c*x))/c)
 genexpon = genexpon_gen(a=0.0,name='genexpon',
                         longname='A generalized exponential',
                         shapes='a,b,c',extradoc="""
 
-Generalized exponential distribution
+Generalized exponential distribution (Ryu 1993)
 
-genexpon.pdf(x,a,b,c) = (a+b*(1-exp(-c*x))) * exp(a*x-b*x+b/c*(1-exp(-c*x)))
+f(x,a,b,c) = (a+b*(1-exp(-c*x))) * exp(-a*x-b*x+b/c*(1-exp(-c*x)))
 for x >= 0, a,b,c > 0.
+
+a, b, c are the first, second and third shape parameters.
+
+References
+----------
+"The Exponential Distribution: Theory, Methods and Applications",
+N. Balakrishnan, Asit P. Basu
 """
                         )
 

Modified: trunk/scipy/stats/tests/test_distributions.py
===================================================================
--- trunk/scipy/stats/tests/test_distributions.py	2008-11-03 19:51:35 UTC (rev 4971)
+++ trunk/scipy/stats/tests/test_distributions.py	2008-11-03 20:38:44 UTC (rev 4972)
@@ -239,6 +239,19 @@
     def test_zero(self):
         assert_equal(stats.expon.pdf(0),1)
 
+class TestGenExpon(TestCase):
+    def test_pdf_unity_area(self):
+        from scipy.integrate import simps
+        # PDF should integrate to one
+        assert_almost_equal(simps(stats.genexpon.pdf(numpy.arange(0,10,0.01),
+                                                     0.5, 0.5, 2.0),
+                                  dx=0.01), 1, 1)
+
+    def test_cdf_bounds(self):
+        # 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 TestDocstring(TestCase):
     def test_docstrings(self):
         """See ticket #761"""