diff options
| -rw-r--r-- | numpy/random/mtrand/distributions.c | 15 | ||||
| -rw-r--r-- | numpy/random/tests/test_regression.py | 8 |
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() |