Allow many distributions to have a scale of 0.

(in which case a stream of 0's is usually returned (or 1's)).

See #5818.

2 files changed, 78 insertions, 60 deletions

diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx
index 3adeb1990..968836ce6 100644
--- a/numpy/random/mtrand/mtrand.pyx
+++ b/numpy/random/mtrand/mtrand.pyx
@@ -1902,15 +1902,13 @@ cdef class RandomState:
         if oloc.shape == oscale.shape == ():
             floc = PyFloat_AsDouble(loc)
             fscale = PyFloat_AsDouble(scale)
-
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
-
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont2_array_sc(self.internal_state, rk_normal, size, floc,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oscale, 0)):
-            raise ValueError("scale <= 0")
+        if np.any(np.less(oscale, 0)):
+            raise ValueError("scale < 0")
         return cont2_array(self.internal_state, rk_normal, size, oloc, oscale,
                            self.lock)
 
@@ -2028,14 +2026,13 @@ cdef class RandomState:
 
         if oscale.shape == ():
             fscale = PyFloat_AsDouble(scale)
-
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont1_array_sc(self.internal_state, rk_exponential, size,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oscale, 0.0)):
-            raise ValueError("scale <= 0")
+        if np.any(np.less(oscale, 0.0)):
+            raise ValueError("scale < 0")
         return cont1_array(self.internal_state, rk_exponential, size, oscale,
                            self.lock)
 
@@ -2146,14 +2143,13 @@ cdef class RandomState:
 
         if oshape.shape == ():
             fshape = PyFloat_AsDouble(shape)
-
-            if fshape <= 0:
-                raise ValueError("shape <= 0")
+            if fshape < 0:
+                raise ValueError("shape < 0")
             return cont1_array_sc(self.internal_state, rk_standard_gamma,
                                   size, fshape, self.lock)
 
-        if np.any(np.less_equal(oshape, 0.0)):
-            raise ValueError("shape <= 0")
+        if np.any(np.less(oshape, 0.0)):
+            raise ValueError("shape < 0")
         return cont1_array(self.internal_state, rk_standard_gamma, size,
                            oshape, self.lock)
 
@@ -2239,18 +2235,17 @@ cdef class RandomState:
         if oshape.shape == oscale.shape == ():
             fshape = PyFloat_AsDouble(shape)
             fscale = PyFloat_AsDouble(scale)
-
-            if fshape <= 0:
-                raise ValueError("shape <= 0")
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
+            if fshape < 0:
+                raise ValueError("shape < 0")
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont2_array_sc(self.internal_state, rk_gamma, size, fshape,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oshape, 0.0)):
-            raise ValueError("shape <= 0")
-        if np.any(np.less_equal(oscale, 0.0)):
-            raise ValueError("scale <= 0")
+        if np.any(np.less(oshape, 0.0)):
+            raise ValueError("shape < 0")
+        if np.any(np.less(oscale, 0.0)):
+            raise ValueError("scale < 0")
         return cont2_array(self.internal_state, rk_gamma, size, oshape, oscale,
                            self.lock)
 
@@ -3121,14 +3116,13 @@ cdef class RandomState:
 
         if oa.shape == ():
             fa = PyFloat_AsDouble(a)
-
-            if fa <= 0:
-                raise ValueError("a <= 0")
+            if fa < 0:
+                raise ValueError("a < 0")
             return cont1_array_sc(self.internal_state, rk_weibull, size, fa,
                                   self.lock)
 
-        if np.any(np.less_equal(oa, 0.0)):
-            raise ValueError("a <= 0")
+        if np.any(np.less(oa, 0.0)):
+            raise ValueError("a < 0")
         return cont1_array(self.internal_state, rk_weibull, size, oa,
                            self.lock)
 
@@ -3234,14 +3228,13 @@ cdef class RandomState:
 
         if oa.shape == ():
             fa = PyFloat_AsDouble(a)
-
-            if fa <= 0:
-                raise ValueError("a <= 0")
+            if fa < 0:
+                raise ValueError("a < 0")
             return cont1_array_sc(self.internal_state, rk_power, size, fa,
                                   self.lock)
 
-        if np.any(np.less_equal(oa, 0.0)):
-            raise ValueError("a <= 0")
+        if np.any(np.less(oa, 0.0)):
+            raise ValueError("a < 0")
         return cont1_array(self.internal_state, rk_power, size, oa, self.lock)
 
     def laplace(self, loc=0.0, scale=1.0, size=None):
@@ -3332,14 +3325,13 @@ cdef class RandomState:
         if oloc.shape == oscale.shape == ():
             floc = PyFloat_AsDouble(loc)
             fscale = PyFloat_AsDouble(scale)
-
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont2_array_sc(self.internal_state, rk_laplace, size, floc,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oscale, 0.0)):
-            raise ValueError("scale <= 0")
+        if np.any(np.less(oscale, 0.0)):
+            raise ValueError("scale < 0")
         return cont2_array(self.internal_state, rk_laplace, size, oloc, oscale,
                            self.lock)
 
@@ -3464,14 +3456,13 @@ cdef class RandomState:
         if oloc.shape == oscale.shape == ():
             floc = PyFloat_AsDouble(loc)
             fscale = PyFloat_AsDouble(scale)
-
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont2_array_sc(self.internal_state, rk_gumbel, size, floc,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oscale, 0.0)):
-            raise ValueError("scale <= 0")
+        if np.any(np.less(oscale, 0.0)):
+            raise ValueError("scale < 0")
         return cont2_array(self.internal_state, rk_gumbel, size, oloc, oscale,
                            self.lock)
 
@@ -3558,14 +3549,13 @@ cdef class RandomState:
         if oloc.shape == oscale.shape == ():
             floc = PyFloat_AsDouble(loc)
             fscale = PyFloat_AsDouble(scale)
-
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont2_array_sc(self.internal_state, rk_logistic, size, floc,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oscale, 0.0)):
-            raise ValueError("scale <= 0")
+        if np.any(np.less(oscale, 0.0)):
+            raise ValueError("scale < 0")
         return cont2_array(self.internal_state, rk_logistic, size, oloc,
                            oscale, self.lock)
 
@@ -3683,14 +3673,13 @@ cdef class RandomState:
         if omean.shape == osigma.shape == ():
             fmean = PyFloat_AsDouble(mean)
             fsigma = PyFloat_AsDouble(sigma)
-
-            if fsigma <= 0:
-                raise ValueError("sigma <= 0")
+            if fsigma < 0:
+                raise ValueError("sigma < 0")
             return cont2_array_sc(self.internal_state, rk_lognormal, size,
                                   fmean, fsigma, self.lock)
 
-        if np.any(np.less_equal(osigma, 0.0)):
-            raise ValueError("sigma <= 0.0")
+        if np.any(np.less(osigma, 0.0)):
+            raise ValueError("sigma < 0.0")
         return cont2_array(self.internal_state, rk_lognormal, size, omean,
                            osigma, self.lock)
 
@@ -3763,14 +3752,13 @@ cdef class RandomState:
 
         if oscale.shape == ():
             fscale = PyFloat_AsDouble(scale)
-
-            if fscale <= 0:
-                raise ValueError("scale <= 0")
+            if fscale < 0:
+                raise ValueError("scale < 0")
             return cont1_array_sc(self.internal_state, rk_rayleigh, size,
                                   fscale, self.lock)
 
-        if np.any(np.less_equal(oscale, 0.0)):
-            raise ValueError("scale <= 0.0")
+        if np.any(np.less(oscale, 0.0)):
+            raise ValueError("scale < 0.0")
         return cont1_array(self.internal_state, rk_rayleigh, size, oscale,
                            self.lock)
 
diff --git a/numpy/random/tests/test_random.py b/numpy/random/tests/test_random.py
index a07fa52f5..3f4952d20 100644
--- a/numpy/random/tests/test_random.py
+++ b/numpy/random/tests/test_random.py
@@ -481,6 +481,9 @@ class TestRandomDist(TestCase):
                             [0.68717433461363442, 1.69175666993575979]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_exponential_0(self):
+        assert_equal(np.random.exponential(scale=0), 0)
+
     def test_f(self):
         np.random.seed(self.seed)
         actual = np.random.f(12, 77, size=(3, 2))
@@ -497,6 +500,9 @@ class TestRandomDist(TestCase):
                             [31.71863275789960568, 33.30143302795922011]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_gamma_0(self):
+        assert_equal(np.random.gamma(shape=0, scale=0), 0)
+
     def test_geometric(self):
         np.random.seed(self.seed)
         actual = np.random.geometric(.123456789, size=(3, 2))
@@ -513,6 +519,9 @@ class TestRandomDist(TestCase):
                             [1.10651090478803416, -0.69535848626236174]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_gumbel_0(self):
+        assert_equal(np.random.gumbel(scale=0), 0)
+
     def test_hypergeometric(self):
         np.random.seed(self.seed)
         actual = np.random.hypergeometric(10.1, 5.5, 14, size=(3, 2))
@@ -547,6 +556,9 @@ class TestRandomDist(TestCase):
                             [-0.05391065675859356, 1.74901336242837324]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_laplace_0(self):
+        assert_equal(np.random.laplace(scale=0), 0)
+
     def test_logistic(self):
         np.random.seed(self.seed)
         actual = np.random.logistic(loc=.123456789, scale=2.0, size=(3, 2))
@@ -555,6 +567,9 @@ class TestRandomDist(TestCase):
                             [-0.21682183359214885, 2.63373365386060332]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_laplace_0(self):
+        assert_(np.random.laplace(scale=0) in [0, 1])
+
     def test_lognormal(self):
         np.random.seed(self.seed)
         actual = np.random.lognormal(mean=.123456789, sigma=2.0, size=(3, 2))
@@ -563,6 +578,9 @@ class TestRandomDist(TestCase):
                             [65.72798501792723869, 86.84341601437161273]])
         assert_array_almost_equal(actual, desired, decimal=13)
 
+    def test_lognormal_0(self):
+        assert_equal(np.random.lognormal(sigma=0), 1)
+
     def test_logseries(self):
         np.random.seed(self.seed)
         actual = np.random.logseries(p=.923456789, size=(3, 2))
@@ -653,6 +671,9 @@ class TestRandomDist(TestCase):
                             [4.18552478636557357, 4.46410668111310471]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_normal_0(self):
+        assert_equal(np.random.normal(scale=0), 0)
+
     def test_pareto(self):
         np.random.seed(self.seed)
         actual = np.random.pareto(a=.123456789, size=(3, 2))
@@ -700,6 +721,9 @@ class TestRandomDist(TestCase):
                             [11.06066537006854311, 17.35468505778271009]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_rayleigh_0(self):
+        assert_equal(np.random.rayleigh(scale=0), 0)
+
     def test_standard_cauchy(self):
         np.random.seed(self.seed)
         actual = np.random.standard_cauchy(size=(3, 2))
@@ -724,6 +748,9 @@ class TestRandomDist(TestCase):
                             [7.54838614231317084, 8.012756093271868]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_standard_gamma_0(self):
+        assert_equal(np.random.standard_gamma(shape=0), 0)
+
     def test_standard_normal(self):
         np.random.seed(self.seed)
         actual = np.random.standard_normal(size=(3, 2))
@@ -799,6 +826,9 @@ class TestRandomDist(TestCase):
                             [0.67057783752390987, 1.39494046635066793]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_weibull_0(self):
+        assert_equal(np.random.weibull(a=0), 0)
+
     def test_zipf(self):
         np.random.seed(self.seed)
         actual = np.random.zipf(a=1.23, size=(3, 2))