BUG: np.random.beta with small parameters produces NaNs

Fixes #5851

2 files changed, 22 insertions, 1 deletions

diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c
index ff936fdd8..f5ee6d8c1 100644
--- a/numpy/random/mtrand/distributions.c
+++ b/numpy/random/mtrand/distributions.c
@@ -199,7 +199,20 @@ double rk_beta(rk_state *state, double a, double b)
 
             if ((X + Y) <= 1.0)
             {
-                return X / (X + Y);
+                if (X +Y > 0)
+                {
+                    return X / (X + Y);
+                }
+                else
+                {
+                    double logX = log(U) / a;
+                    double logY = log(V) / b;
+                    double logM = logX > logY ? logX : logY;
+                    logX -= logM;
+                    logY -= logM;
+
+                    return exp(logX - log(exp(logX) + exp(logY)));
+                }
             }
         }
     }
diff --git a/numpy/random/tests/test_regression.py b/numpy/random/tests/test_regression.py
index 1a5854e82..52be0d4d5 100644
--- a/numpy/random/tests/test_regression.py
+++ b/numpy/random/tests/test_regression.py
@@ -93,5 +93,13 @@ class TestRegression(TestCase):
         np.random.multivariate_normal([0], [[0]], size=np.int_(1))
         np.random.multivariate_normal([0], [[0]], size=np.int64(1))
 
+    def test_beta_small_parameters(self):
+        # Test that beta with small a and b parameters does not produce
+        # NaNs due to roundoff errors causing 0 / 0, gh-5851
+        np.random.seed(1234567890)
+        x = np.random.beta(0.0001, 0.0001, size=100)
+        assert_(not np.any(np.isnan(x)), 'Nans in np.random.beta')
+
+
 if __name__ == "__main__":
     run_module_suite()