BUG: Fixes random.noncentral_chisquare when 0 < df <= 1

Closes #5766.

3 files changed, 21 insertions, 8 deletions

diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c
index 84174e105..ff936fdd8 100644
--- a/numpy/random/mtrand/distributions.c
+++ b/numpy/random/mtrand/distributions.c
@@ -218,11 +218,17 @@ double rk_chisquare(rk_state *state, double df)
 
 double rk_noncentral_chisquare(rk_state *state, double df, double nonc)
 {
-    double Chi2, N;
-
-    Chi2 = rk_chisquare(state, df-1);
-    N = rk_gauss(state) + sqrt(nonc);
-    return Chi2 + N*N;
+    if(1 < df)
+    {
+        const double Chi2 = rk_chisquare(state, df - 1);
+        const double N = rk_gauss(state) + sqrt(nonc);
+        return Chi2 + N*N;
+    }
+    else
+    {
+        const long i = rk_poisson(state, nonc / 2.0);
+        return rk_chisquare(state, df + 2 * i);
+    }
 }
 
 double rk_f(rk_state *state, double dfnum, double dfden)
diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx
index c4927a3f3..69fc0414b 100644
--- a/numpy/random/mtrand/mtrand.pyx
+++ b/numpy/random/mtrand/mtrand.pyx
@@ -2201,7 +2201,8 @@ cdef class RandomState:
         Parameters
         ----------
         df : int
-            Degrees of freedom, should be >= 1.
+            Degrees of freedom, should be > 0 as of Numpy 1.10,
+            should be > 1 for earlier versions.
         nonc : float
             Non-centrality, should be > 0.
         size : int or tuple of ints, optional
@@ -2267,7 +2268,7 @@ cdef class RandomState:
         fdf = PyFloat_AsDouble(df)
         fnonc = PyFloat_AsDouble(nonc)
         if not PyErr_Occurred():
-            if fdf <= 1:
+            if fdf <= 0:
                 raise ValueError("df <= 0")
             if fnonc <= 0:
                 raise ValueError("nonc <= 0")
@@ -2279,7 +2280,7 @@ cdef class RandomState:
         odf = <ndarray>PyArray_FROM_OTF(df, NPY_DOUBLE, NPY_ARRAY_ALIGNED)
         ononc = <ndarray>PyArray_FROM_OTF(nonc, NPY_DOUBLE, NPY_ARRAY_ALIGNED)
         if np.any(np.less_equal(odf, 0.0)):
-            raise ValueError("df <= 1")
+            raise ValueError("df <= 0")
         if np.any(np.less_equal(ononc, 0.0)):
             raise ValueError("nonc < 0")
         return cont2_array(self.internal_state, rk_noncentral_chisquare, size,
diff --git a/numpy/random/tests/test_random.py b/numpy/random/tests/test_random.py
index 897a8e6f0..596c218a2 100644
--- a/numpy/random/tests/test_random.py
+++ b/numpy/random/tests/test_random.py
@@ -498,6 +498,12 @@ class TestRandomDist(TestCase):
                             [5.03461598262724586, 17.94973089023519464]])
         np.testing.assert_array_almost_equal(actual, desired, decimal=14)
 
+        actual = np.random.noncentral_chisquare(df=.5, nonc=.2, size=(3, 2))
+        desired = np.array([[ 1.47145377828516666,  0.15052899268012659],
+                            [ 0.00943803056963588,  1.02647251615666169],
+                            [ 0.332334982684171  ,  0.15451287602753125]])
+        np.testing.assert_array_almost_equal(actual, desired, decimal=14)
+
     def test_noncentral_f(self):
         np.random.seed(self.seed)
         actual = np.random.noncentral_f(dfnum=5, dfden=2, nonc=1,